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_MatrixVarReductImplicit.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_MatrixVarReductImplicit.hpp"
43 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
44 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
45 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
46 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
47 #include "AbstractLinAlgPack_SpVectorClass.hpp"
48 #include "AbstractLinAlgPack_AssertOp.hpp"
49 #include "Teuchos_Workspace.hpp"
50 #include "Teuchos_Assert.hpp"
51 
52 namespace {
53 
54 //
55 // Implicit matrix-vector multiplication:
56 //
57 // y = b*y + a*op(inv(C)*N)*x
58 //
59 template<class V>
60 void imp_Vp_StMtV_implicit(
62  ,AbstractLinAlgPack::value_type a
65  ,BLAS_Cpp::Transp D_trans
66  ,const V &x
67  ,DenseLinAlgPack::value_type b
68  )
69 {
70  using BLAS_Cpp::no_trans;
71  using BLAS_Cpp::trans;
72  namespace alap = AbstractLinAlgPack;
73  using Teuchos::Workspace;
75 
77  r = C.rows(),
78  dof = N.cols();
79 
80  // ToDo: Pass in workspace vectors to save some allocation time!
81 
82  if( D_trans == no_trans ) {
83  // y = b*y
84  alap::Vt_S(y,b);
85  //
86  // y += -a * inv(C) * ( N * x )
87  //
88  alap::VectorSpace::vec_mut_ptr_t
89  t1 = N.space_cols().create_member(),
90  t2 = C.space_rows().create_member();
91  // t1 = N*x
92  LinAlgOpPack::V_MtV( t1.get(), N, no_trans, x );
93  // t2 = inv(C) * t1
94  alap::V_InvMtV( t2.get(), C, no_trans, *t1 );
95  // y += a*t2
96  alap::Vp_StV( y, -a, *t2 );
97  }
98  else {
99  //
100  // y = b*y - a * N' * ( inv(C') * x )
101  //
102  alap::VectorSpace::vec_mut_ptr_t
103  t1 = C.space_cols().create_member();
104  // t1 = inv(C')*x
105  alap::V_InvMtV( t1.get(), C, trans, x );
106  // y = b*y - a*N'*t1
107  alap::Vp_StMtV( y, -a, N, trans, *t1, b );
108  }
109 }
110 
111 /*
112 
113 //
114 // Generate a row of inv(C)*N if not already computed.
115 //
116 void imp_compute_InvCtN_row(
117  DenseLinAlgPack::size_type r
118  ,const ConstrainedOptPack::DecompositionSystemVarReduct &decomp_sys
119  ,DenseLinAlgPack::size_type j
120  ,DenseLinAlgPack::DVectorSlice *e_j // Set to all zeros on input and output!
121  ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t *InvCtN_rows
122  )
123 {
124  typedef DenseLinAlgPack::value_type value_type;
125  using DenseLinAlgPack::DVectorSlice;
126  if( (*InvCtN_rows)[j-1] == NULL ) {
127  // Generate row j of inv(C)*N
128  value_type *vec = (*InvCtN_rows)[j-1] = new value_type[r]; // ToDo: We may want to allocate more vectors at once!
129  DVectorSlice row_j(vec,r);
130  // row_j = N'*inv(C')*e_j
131  (*e_j)(j) = 1.0;
132  imp_Vp_StMtV_implicit( &row_j, 1.0, decomp_sys, BLAS_Cpp::trans, *e_j, 0.0 );
133  (*e_j)(j) = 0.0;
134  }
135 }
136 
137 //
138 // Perform the matrix-vector multiplication:
139 //
140 // y = b*y -a * op(P) * [inv(C) * N] * x
141 //
142 // by generating rows [inv(C)*N](j,:) for each nonzero entry op(P)(i,j).
143 //
144 template<class V>
145 void imp_Vp_StPtMtV_by_row(
146  DenseLinAlgPack::DVectorSlice *y
147  ,DenseLinAlgPack::value_type a
148  ,const ConstrainedOptPack::GenPermMatrixSlice &P
149  ,BLAS_Cpp::Transp P_trans
150  ,const ConstrainedOptPack::DecompositionSystemVarReduct &decomp_sys
151  ,const V &x
152  ,DenseLinAlgPack::value_type b
153  ,ConstrainedOptPack::MatrixVarReductImplicit::InvCtN_rows_t *InvCtN_rows
154  )
155 {
156  using BLAS_Cpp::no_trans;
157  using BLAS_Cpp::trans;
158  using BLAS_Cpp::trans_not;
159  using BLAS_Cpp::rows;
160  using BLAS_Cpp::cols;
161  using DenseLinAlgPack::dot;
162  using DenseLinAlgPack::DVectorSlice;
163  using AbstractLinAlgPack::dot;
164  using AbstractLinAlgPack::Vp_StMtV;
165  using AbstractLinAlgPack::GenPermMatrixSlice;
166  using Teuchos::Workspace;
167  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
168 
169  const DenseLinAlgPack::size_type
170  D_rows = decomp_sys.C().rows(),
171  D_cols = decomp_sys.N().cols();
172  // y = b*y
173  if(b==0.0) *y = 0.0;
174  else if(b!=1.0) DenseLinAlgPack::Vt_S(y,b);
175  // Compute t = N'*inv(C')*e(j) then y(i) += -a*t'*x where op(P)(i,j) = 1.0
176  Workspace<DenseLinAlgPack::value_type> e_j_ws(wss,D_rows);
177  DVectorSlice e_j(&e_j_ws[0],e_j_ws.size());
178  e_j = 0.0;
179  for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
180  const DenseLinAlgPack::size_type
181  i = P_trans == no_trans ? itr->row_i() : itr->col_j(),
182  j = P_trans == no_trans ? itr->col_j() : itr->row_i();
183  // t = op(M') * e(j)
184  imp_compute_InvCtN_row(D_rows,decomp_sys,j,&e_j,InvCtN_rows);
185  DVectorSlice t((*InvCtN_rows)[j-1],D_cols);
186  // y(i) += -a*t'*x
187  (*y)(i) += (-a) * dot(t,x);
188  }
189 }
190 
191 */
192 
193 } // end namespace
194 
195 namespace ConstrainedOptPack {
196 
198  const mat_nonsing_ptr_t &C
199  ,const mat_ptr_t &N
200  ,const mat_ptr_t &D_direct
201  )
202 {
203  namespace rcp = MemMngPack;
204  // Validate the inputs
206  C.get() == NULL, std::invalid_argument
207  ,"MatrixVarReductImplicit::initialize(...): Error, "
208  "C.get() must not be NULL" );
210  N.get() == NULL, std::invalid_argument
211  ,"MatrixVarReductImplicit::initialize(...): Error, "
212  "N.get() must not be NULL" );
213  if( D_direct.get() ) {
214  const bool is_compatible_cols = D_direct->space_cols().is_compatible(C->space_cols());
216  !is_compatible_cols, VectorSpace::IncompatibleVectorSpaces
217  ,"MatrixVarReductImplicit::initialize(...): Error, "
218  "D_direct->space_cols() is not compatible with C->space_cols()" );
219  const bool is_compatible_rows = D_direct->space_rows().is_compatible(N->space_rows());
221  !is_compatible_rows, VectorSpace::IncompatibleVectorSpaces
222  ,"MatrixVarReductImplicit::initialize(...): Error, "
223  "D_direct->space_rows() is not compatible with N->space_rows()" );
224  }
225  // Set the members
226  C_ = C;
227  N_ = N;
228  D_direct_ = D_direct;
229  if(!InvCtN_rows_set_list_.empty()) { // Free previously allocated vectors
230  for( InvCtN_rows_set_list_t::iterator itr = InvCtN_rows_set_list_.begin();
231  itr != InvCtN_rows_set_list_.end(); ++itr )
232  {
233  InvCtN_rows_[*itr] = Teuchos::null;
234  }
235  InvCtN_rows_set_list_.clear();
236  }
237 }
238 
240 {
241  namespace rcp = MemMngPack;
242  C_ = Teuchos::null;
243  N_ = Teuchos::null;
244  D_direct_ = Teuchos::null;
245 }
246 
247 // Overridden from MatrixBase
248 
250 {
251  return C_.get() ? C_->rows() : 0;
252 }
253 
255 {
256  return N_.get() ? N_->cols() : 0;
257 }
258 
259 // Overridden from MatrixOp
260 
261 const VectorSpace& MatrixVarReductImplicit::space_cols() const
262 {
263  assert_initialized();
264  return C_->space_cols();
265 }
266 
267 const VectorSpace& MatrixVarReductImplicit::space_rows() const
268 {
269  assert_initialized();
270  return N_->space_rows();
271 }
272 
273 MatrixOp& MatrixVarReductImplicit::operator=(const MatrixOp& M)
274 {
275  assert_initialized();
276  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
277  return *this;
278 }
279 
280 std::ostream& MatrixVarReductImplicit::output(std::ostream& o) const
281 {
282  o << "This is a " << this->rows() << " x " << this->cols()
283  << " variable reduction matrix D = -inv(C)*N where C and N are:\n"
284  << "C =\n" << *C_
285  << "N =\n" << *N_;
286  return o;
287 }
288 
290  VectorMutable* y, value_type a
291  ,BLAS_Cpp::Transp D_trans
292  ,const Vector& x, value_type b
293  ) const
294 {
295  assert_initialized();
297  imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
298 }
299 
301  VectorMutable* y, value_type a
302  ,BLAS_Cpp::Transp D_trans
303  ,const SpVectorSlice& x, value_type b
304  ) const
305 {
306  using BLAS_Cpp::rows;
307  using BLAS_Cpp::cols;
308  using Teuchos::Workspace;
310 
311  assert_initialized();
313  imp_Vp_StMtV_implicit( y, a, *C_, *N_, D_trans, x, b );
314 /*
315  const size_type
316  D_rows = this->rows(), D_cols = this->cols(),
317  opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
318  DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),D_rows,D_cols,D_trans,x.size());
319  if( use_dense_mat_vec_ && D_dense_.rows() > 0 ) {
320  LinAlgOpPack::Vp_StMtV( y, a, D_dense_, D_trans, x, b );
321  }
322  else {
323  if( x.nz() == x.size() ) { // This is B.S. Should use MatrixWithOpFactorized object for C!
324  DVectorSlice dx = AbstractLinAlgPack::dense_view(x);
325  imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx, b );
326  }
327  else if( D_trans == BLAS_Cpp::trans && x.nz() < D_cols ) {
328  //
329  // y = b*y + (-a)*[N'*inv(C')]*x
330  //
331  // We can do something crafty here. We can generate columns of N'*inv(C')
332  // and then perform y += -a*[N'*inv(C')](:,j)*x(j) for nonzero x(j)
333  //
334  Workspace<DenseLinAlgPack::value_type> e_j_ws(wss,D_rows);
335  DVectorSlice e_j(&e_j_ws[0],e_j_ws.size());
336  e_j = 0.0;
337  // y = b*y
338  if(b==0.0) *y = 0.0;
339  else if(b!=1.0) Vt_S(y,b);
340  // y += -a*[N'*inv(C')](:,j)*x(j), for j <: { j | x(j) != 0.0 }
341  const SpVectorSlice::difference_type o = x.offset();
342  for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
343  const size_type j = itr->indice() + o;
344  imp_compute_InvCtN_row(D_rows,*decomp_sys_,j,&e_j,&InvCtN_rows_);
345  DenseLinAlgPack::Vp_StV( y, -a * itr->value(), DVectorSlice(InvCtN_rows_[j-1],D_cols) );
346  }
347  }
348  else { // This is B.S. Should use MatrixWithOpFactorized object for C!
349  DVector dx;
350  LinAlgOpPack::assign( &dx, x );
351  imp_Vp_StMtV_implicit( y, -a, *decomp_sys_, D_trans, dx(), b );
352  }
353  }
354 */
355  // ToDo: Consider implementing the above specialized implementation!
356 }
357 
359  VectorMutable* y, value_type a
360  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
361  ,BLAS_Cpp::Transp D_trans
362  ,const Vector& x, value_type b
363  ) const
364 {
365  using BLAS_Cpp::rows;
366  using BLAS_Cpp::cols;
367  assert_initialized();
368 /*
369  const size_type
370  D_rows = this->rows(), D_cols = this->cols(),
371  opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
372  DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
373  DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
374  if( D_dense_.rows() > 0 ) {
375  AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
376  }
377  else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
378  // Just use the default implementation
379  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
380  }
381  else {
382  imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
383  }
384 */
385  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
386 }
387 
389  VectorMutable* y, value_type a
390  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
391  ,BLAS_Cpp::Transp D_trans
392  ,const SpVectorSlice& x, value_type b
393  ) const
394 {
395  using BLAS_Cpp::rows;
396  using BLAS_Cpp::cols;
397  assert_initialized();
398 /*
399  const size_type
400  D_rows = this->rows(), D_cols = this->cols(),
401  opD_rows = rows(D_rows,D_cols,D_trans), opD_cols = cols(D_rows,D_cols,D_trans);
402  DenseLinAlgPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans,opD_rows);
403  DenseLinAlgPack::Vp_MtV_assert_sizes(cols(P.rows(),P.cols(),P_trans),D_rows,D_cols,D_trans,x.size());
404  if( D_dense_.rows() > 0 ) {
405  AbstractLinAlgPack::dense_Vp_StPtMtV(y,a,P,P_trans,D_dense_,D_trans,x,b);
406  }
407  else if( P.nz() > D_cols || D_trans == BLAS_Cpp::trans ) {
408  // Just use the default implementation
409  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b);
410  }
411  else {
412  imp_Vp_StPtMtV_by_row(y,a,P,P_trans,*decomp_sys_,x,b,&InvCtN_rows_);
413  }
414 */
415  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,D_trans,x,b); // ToDo:Update specialized implementation above!
416 }
417 
418 // Private member functions
419 
420 void MatrixVarReductImplicit::assert_initialized() const
421 {
423  C_.get() == NULL, std::logic_error
424  ,"MatrixVarReductImplicit::assert_initialized(): Error, "
425  "initialize(...) has not been called yet!" );
426 }
427 
428 } // end namespace ConstrainedOptPack
virtual const VectorSpace & space_rows() const =0
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
virtual void initialize(const mat_nonsing_ptr_t &C, const mat_ptr_t &N, const mat_ptr_t &D_direct)
Initialize this matrix object.
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
T * get() const
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta) const
virtual size_type cols() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
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)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
virtual const VectorSpace & space_cols() const =0
void Vp_MtV_assert_compatibility(VectorMutable *v_lhs, const MatrixOp &m_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
virtual vec_mut_ptr_t create_member() const =0
virtual size_type rows() const
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
virtual void set_uninitialized()
Set the matrix to uninitialized.