MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_BasisSystemComposite.cpp
Go to the documentation of this file.
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 <assert.h>
43 
52 #include "Teuchos_dyn_cast.hpp"
53 #include "Teuchos_Assert.hpp"
54 
55 namespace {
56 
57 // Allocates a MultiVectorMutable object given a vector space
58 // object and the number of columns (num_vecs).
59 
60 class AllocatorMultiVectorMutable {
61 public:
62  AllocatorMultiVectorMutable(
65  )
66  :vec_space_(vec_space)
67  ,num_vecs_(num_vecs)
68  {}
69  typedef Teuchos::RCP<
71  ptr_t allocate() const
72  {
73  return vec_space_->create_members(num_vecs_);
74  }
75 private:
78 }; // end class AllocatorMultiVectorMutable
79 
80 } // end namespace
81 
82 namespace AbstractLinAlgPack {
83 
84 // Static member functions
85 
87  const VectorSpace::space_ptr_t &space_xD
88  ,const VectorSpace::space_ptr_t &space_xI
89  ,Range1D *var_dep
90  ,Range1D *var_indep
91  ,VectorSpace::space_ptr_t *space_x
92  )
93 {
94  namespace mmp = MemMngPack;
95 #ifdef TEUCHOS_DEBUG
97  space_xD.get() == NULL, std::invalid_argument
98  ,"BasisSystemComposite::initialize_space_x(...): Error!" );
100  var_dep == NULL, std::invalid_argument
101  ,"BasisSystemComposite::initialize_space_x(...): Error!" );
102  TEUCHOS_TEST_FOR_EXCEPTION(
103  space_xI.get() != NULL && var_indep == NULL, std::invalid_argument
104  ,"BasisSystemComposite::initialize_space_x(...): Error!" );
105  TEUCHOS_TEST_FOR_EXCEPTION(
106  space_x == NULL, std::invalid_argument
107  ,"BasisSystemComposite::initialize_space_x(...): Error!" );
108 #endif
109  *var_dep = Range1D(1,space_xD->dim());
110  if(space_xI.get()) *var_indep = Range1D(var_dep->ubound()+1,var_dep->ubound()+space_xI->dim());
111  else *var_indep = Range1D::Invalid;
112  if(space_xI.get()) {
114  vec_spaces[2] = { space_xD, space_xI };
115  *space_x = Teuchos::rcp(new VectorSpaceBlocked(vec_spaces,2));
116  }
117  else {
118  *space_x = space_xD;
119  }
120 }
121 
124 {
125  namespace mmp = MemMngPack;
127 }
128 
130  const VectorSpace::space_ptr_t &space_x
131  ,const Range1D &var_dep
132  ,const Range1D &var_indep
133  ,const VectorSpace::space_ptr_t &space_c
134  ,const C_ptr_t &C
135  ,const N_ptr_t &N
136  ,MatrixOp *Gc
137  )
138 {
139  namespace mmp = MemMngPack;
140  using Teuchos::dyn_cast;
141 #ifdef TEUCHOS_DEBUG
143  space_x.get() == NULL, std::invalid_argument
144  ,"BasisSystemComposite::initialize_Gc(...): Error!" );
146  space_c.get() == NULL, std::invalid_argument
147  ,"BasisSystemComposite::initialize_Gc(...): Error!" );
148 #endif
149  const size_type
150  n = space_x->dim(),
151  m = space_c->dim(),
152  var_dep_size = var_dep.size();
153 #ifdef TEUCHOS_DEBUG
155  C.get() == NULL, std::invalid_argument
156  ,"BasisSystemComposite::initialize_Gc(...): Error!" );
158  var_dep_size < n && N.get() == NULL, std::invalid_argument
159  ,"BasisSystemComposite::initialize_Gc(...): Error!" );
161  Gc == NULL, std::invalid_argument
162  ,"BasisSystemComposite::initialize_Gc(...): Error!" );
163 #endif
164 
166  &Gc_comp = dyn_cast<MatrixComposite>(*Gc);
167 
168  //
169  // Gc = [ C'; N' ]
170  //
171 
172  Gc_comp.reinitialize(n,m);
173  // Add the C matrix object
175  C_rr_ptr_ptr_t
176  C_rr_ptr_ptr = Teuchos::rcp(new mmp::ReleaseResource_ref_count_ptr<MatrixOpNonsing>(C));
177  Gc_comp.add_matrix(
178  var_dep.lbound()-1, 0 // row_offset, col_offset
179  ,1.0 // alpha
180  ,C_rr_ptr_ptr->ptr.get() // A
181  ,C_rr_ptr_ptr // A_release
182  ,BLAS_Cpp::trans // A_trans
183  );
184  if( n > m ) {
185  // Add the N matrix object
187  N_rr_ptr_ptr_t
188  N_rr_ptr_ptr = Teuchos::rcp(new mmp::ReleaseResource_ref_count_ptr<MatrixOp>(N));
189  Gc_comp.add_matrix(
190  var_indep.lbound()-1, 0 // row_offset, col_offset
191  ,1.0 // alpha
192  ,N_rr_ptr_ptr->ptr.get() // A
193  ,N_rr_ptr_ptr // A_release
194  ,BLAS_Cpp::trans // A_trans
195  );
196  }
197  // Finish construction
198  Gc_comp.finish_construction( space_x, space_c );
199 }
200 
202  MatrixOp *Gc
203  ,MatrixOpNonsing **C
204  ,MatrixOp **N
205  )
206 {
207  using Teuchos::dyn_cast;
208 #ifdef TEUCHOS_DEBUG
210  Gc == NULL, std::invalid_argument
211  ,"BasisSystemComposite::get_C_N(...): Error!" );
212 #endif
213  const size_type
214  n = Gc->rows(),
215  m = Gc->cols();
216 #ifdef TEUCHOS_DEBUG
218  C == NULL, std::invalid_argument
219  ,"BasisSystemComposite::get_C_N(...): Error!" );
221  n > m && N == NULL, std::invalid_argument
222  ,"BasisSystemComposite::get_C_N(...): Error!" );
223 #endif
224  // Get reference to concrete Gc matrix subclass
226  &Gc_comp = dyn_cast<MatrixComposite>(*Gc);
227  // Get referencs to the aggregate C and N matrtices
228  MatrixComposite::matrix_list_t::const_iterator
229  mat_itr = Gc_comp.matrices_begin(),
230  mat_end = Gc_comp.matrices_end();
231  if( mat_itr != mat_end ) {
232  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
233  *C = &dyn_cast<MatrixOpNonsing>(
234  const_cast<MatrixOp&>(*(mat_itr++)->A_) );
235  if( n > m ) {
236  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
237  *N = &const_cast<MatrixOp&>(*(mat_itr++)->A_);
238  }
239  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr == mat_end ) );
240  }
241  else {
242  *C = NULL;
243  *N = NULL;
244  }
245 }
246 
248  const MatrixOp &Gc
249  ,const MatrixOpNonsing **C
250  ,const MatrixOp **N
251  )
252 {
253  using Teuchos::dyn_cast;
254  const size_type
255  n = Gc.rows(),
256  m = Gc.cols();
257 #ifdef TEUCHOS_DEBUG
259  C == NULL, std::invalid_argument
260  ,"BasisSystemComposite::get_C_N(...): Error!" );
262  n > m && N == NULL, std::invalid_argument
263  ,"BasisSystemComposite::get_C_N(...): Error!" );
264 #endif
265  // Get reference to concrete Gc matrix subclass
267  &Gc_comp = dyn_cast<const AbstractLinAlgPack::MatrixComposite>(Gc);
268  // Get referencs to the aggregate C and N matrtices
269  MatrixComposite::matrix_list_t::const_iterator
270  mat_itr = Gc_comp.matrices_begin(),
271  mat_end = Gc_comp.matrices_end();
272  if( mat_itr != mat_end ) {
273  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
274  *C = &dyn_cast<const MatrixOpNonsing>(*(mat_itr++)->A_);
275  if( n > m ) {
276  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
277  *N = &dyn_cast<const MatrixOp>(*(mat_itr++)->A_);
278  }
279  TEUCHOS_TEST_FOR_EXCEPT( !( mat_itr == mat_end ) );
280  }
281  else {
283  true, std::invalid_argument
284  ,"BasisSystemComposite::get_C_N(...): Error, "
285  "The Gc matrix object has not been initialized with C and N!" );
286  }
287 }
288 
289 // Constructors / initializers
290 
292  :BasisSystem(Teuchos::null,Teuchos::null)
293 {}
294 
296  const VectorSpace::space_ptr_t &space_x
297  ,const VectorSpace::space_ptr_t &space_c
298  ,const mat_nonsing_fcty_ptr_t &factory_C
299  ,const mat_sym_fcty_ptr_t &factory_transDtD
300  ,const mat_sym_nonsing_fcty_ptr_t &factory_S
301  )
302  :BasisSystem(Teuchos::null,Teuchos::null)
303 {
304  namespace mmp = MemMngPack;
305  const size_type n = space_x->dim(), m = space_c->dim();
306  this->initialize(
307  space_x,Range1D(1,space_c->dim()),(n>m ? Range1D(space_c->dim()+1,space_x->dim()) : Range1D::Invalid)
309  );
310 }
311 
313  const VectorSpace::space_ptr_t &space_x
314  ,const Range1D &var_dep
315  ,const Range1D &var_indep
316  ,const VectorSpace::space_ptr_t &space_c
317  ,const mat_nonsing_fcty_ptr_t &factory_C
318  ,const mat_sym_fcty_ptr_t &factory_transDtD
319  ,const mat_sym_nonsing_fcty_ptr_t &factory_S
320  ,const mat_fcty_ptr_t &factory_D
321  )
322  :BasisSystem(Teuchos::null,Teuchos::null)
323 {
324  this->initialize(
325  space_x,var_dep,var_indep,space_c,factory_C,factory_transDtD,factory_S
326  ,factory_D
327  );
328 }
329 
331  const VectorSpace::space_ptr_t &space_x
332  ,const Range1D &var_dep
333  ,const Range1D &var_indep
334  ,const VectorSpace::space_ptr_t &space_c
335  ,const mat_nonsing_fcty_ptr_t &factory_C
336  ,const mat_sym_fcty_ptr_t &factory_transDtD
337  ,const mat_sym_nonsing_fcty_ptr_t &factory_S
338  ,const mat_fcty_ptr_t &factory_D
339  )
340 {
341  namespace mmp = MemMngPack;
342 #ifdef TEUCHOS_DEBUG
344  space_x.get() == NULL, std::invalid_argument
345  ,"BasisSystemComposite::initialize(...): Error!" );
347  space_c.get() == NULL, std::invalid_argument
348  ,"BasisSystemComposite::initialize(...): Error!" );
349 #endif
350  const size_type n = space_x->dim(), m = space_c->dim();
351 #ifdef TEUCHOS_DEBUG
353  var_dep.size() + var_indep.size() != space_x->dim(), std::invalid_argument
354  ,"BasisSystemComposite::initialize(...): Error!" );
356  n > m && var_dep.lbound() < var_indep.lbound() && (var_dep.lbound() != 1 || var_dep.ubound()+1 != var_indep.lbound())
357  , std::invalid_argument
358  ,"BasisSystemComposite::initialize(...): Error!" );
359  TEUCHOS_TEST_FOR_EXCEPTION(
360  n > m && var_dep.lbound() >= var_indep.lbound() && (var_indep.lbound() != 1 || var_indep.ubound()+1 != var_dep.lbound())
361  , std::invalid_argument
362  ,"BasisSystemComposite::initialize(...): Error!" );
363  TEUCHOS_TEST_FOR_EXCEPTION(
364  factory_C.get() == NULL, std::invalid_argument
365  ,"BasisSystemComposite::initialize(...): Error!" );
366 #endif
367 
368  space_x_ = space_x;
369  var_dep_ = var_dep;
370  var_indep_ = var_indep;
371  space_c_ = space_c;
372  factory_C_ = factory_C;
373  factory_D_ = factory_D;
374  if( factory_D_.get() == NULL ) {
375  factory_D_ = Teuchos::abstractFactoryStd<MatrixOp,MultiVectorMutable>(
376  AllocatorMultiVectorMutable(space_x_->sub_space(var_dep),var_indep.size() ) );
377  }
378  BasisSystem::initialize(factory_transDtD,factory_S);
379 }
380 
382 {
383  namespace mmp = MemMngPack;
384 
385  space_x_ = Teuchos::null;
386  var_dep_ = Range1D::Invalid;
387  var_indep_ = Range1D::Invalid;
388  factory_C_ = Teuchos::null;
389  factory_D_ = Teuchos::null;
390 }
391 
394 {
395  return space_x_;
396 }
397 
400 {
401  return space_c_;
402 }
403 
404 // To be overridden by subclasses
405 
407  const MatrixOpNonsing &C
408  ,const MatrixOp &N
409  ,MatrixOp *D
410  ,EMatRelations mat_rel
411  ) const
412 {
414  M_StInvMtM( D, -1.0, C, BLAS_Cpp::no_trans, N, BLAS_Cpp::no_trans ); // D = -inv(C)*N
415 }
416 
417 // Overridden from BasisSytem
418 
421 {
422  return factory_C_;
423 }
424 
427 {
428  return factory_D_;
429 }
430 
432 {
433  return var_dep_;
434 }
435 
437 {
438  return var_indep_;
439 }
440 
442  const MatrixOp &Gc
444  ,MatrixOp *D
445  ,MatrixOp *GcUP
446  ,EMatRelations mat_rel
447  ,std::ostream *out
448  ) const
449 {
450  using Teuchos::dyn_cast;
451  const index_type
452  n = var_dep_.size() + var_indep_.size(),
453  m = var_dep_.size();
455  n == 0, std::logic_error
456  ,"BasisSystemComposite::update_basis(...): Error, this must be initialized first!" );
458  C == NULL && ( n > m ? D == NULL : false ), std::logic_error
459  ,"BasisSystemComposite::update_basis(...): Error, C or D must be non-NULL!" );
460  // Get references to the aggregate C and N matrices
461  const MatrixOpNonsing
462  *C_aggr = NULL;
463  const MatrixOp
464  *N_aggr = NULL;
465  get_C_N( Gc, &C_aggr, &N_aggr );
466  // Setup C
467  if( C ) {
468  *C = *C_aggr; // Assignment had better work!
469  }
470  // Compute D
471  if( D ) {
472  update_D(*C_aggr,*N_aggr,D,mat_rel);
473  }
474 }
475 
476 } // end namespace AbstractLinAlgPack
static void get_C_N(MatrixOp *Gc, MatrixOpNonsing **C, MatrixOp **N)
Get the non-const aggregate matrices C and N (or NULL pointers if not initialized).
Interface for the creation and maintainance of a basis matrix for a decomposition of linearlized cons...
static void initialize_Gc(const VectorSpace::space_ptr_t &space_x, const Range1D &var_dep, const Range1D &var_indep, const VectorSpace::space_ptr_t &space_c, const C_ptr_t &C, const N_ptr_t &N, MatrixOp *Gc)
Initialize the Gc matrix object given created from space_Gc()->create().
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index size() const
Return the size of the range (ubound() - lbound() + 1)
void add_matrix(size_type row_offset, size_type col_offset, value_type alpha, const GenPermMatrixSlice *P, const release_resource_ptr_t &P_release, BLAS_Cpp::Transp P_trans, const MatrixOp *A, const release_resource_ptr_t &A_release, BLAS_Cpp::Transp A_trans, const GenPermMatrixSlice *Q, const release_resource_ptr_t &Q_release, BLAS_Cpp::Transp Q_trans)
Add a sub-matrix alpha*op(P)*op(A)*op(Q).
Index ubound() const
Return upper bound of the range.
virtual void update_D(const MatrixOpNonsing &C, const MatrixOp &N, MatrixOp *D, EMatRelations mat_rel) const
Overridden by subclasses to update D if a specialized implementation is needed.
T * get() const
void update_basis(const MatrixOp &Gc, MatrixOpNonsing *C, MatrixOp *D, MatrixOp *GcUP, EMatRelations mat_rel, std::ostream *out) const
Transposed.
static void initialize_space_x(const VectorSpace::space_ptr_t &space_xD, const VectorSpace::space_ptr_t &space_xI, Range1D *var_dep, Range1D *var_indep, VectorSpace::space_ptr_t *space_x)
Initialize the composite vector space for x = [ xD; xI ] as well as var_dep and var_indep.
virtual void finish_construction(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows)
Call to finish the construction process.
Not transposed.
virtual size_type cols() const
Return the number of columns in the matrix.
virtual void initialize(const mat_sym_fcty_ptr_t &factory_transDtD, const mat_sym_nonsing_fcty_ptr_t &factory_S)
Initialize the factory objects for the special matrices for D'*D and S = I + D'*D.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
. One-based subregion index range class.
T_To & dyn_cast(T_From &from)
std::ostream * out
VectorSpace subclass for the composite of one or more VectorSpace objects.
Base class for all matrices that support basic matrix operations.
static const Range1D Invalid
Range1D(INVALID)
void reinitialize(size_type rows=0, size_type cols=0)
Initialize a sized (on unsized) zero matrix to start with.
Index lbound() const
Return lower bound of the range.
Interface for a collection of mutable vectors (multi-vector, matrix).
virtual const mat_sym_nonsing_fcty_ptr_t factory_S() const
Returns a matrix factory for the result of S = I + D'*D
void initialize(const VectorSpace::space_ptr_t &space_x, const Range1D &var_dep, const Range1D &var_indep, const VectorSpace::space_ptr_t &space_c, const mat_nonsing_fcty_ptr_t &factory_C, const mat_sym_fcty_ptr_t &factory_transDtD, const mat_sym_nonsing_fcty_ptr_t &factory_S, const mat_fcty_ptr_t &factory_D=Teuchos::null)
Initialize.
virtual const mat_sym_fcty_ptr_t factory_transDtD() const
Returns a matrix factory for the result of J = D'*D
Abstract base class for all nonsingular polymorphic matrices that can be used to compute matrix-vecto...
const f_int const f_int & N
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
virtual size_type rows() const
Return the number of rows in the matrix.
Matrix class for matrices composed out of a set of other matrices and vectors.
int n
static const fcty_Gc_ptr_t factory_Gc()
Return a matrix factory object for the composte Gc matrix object.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)