AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_MatrixComposite.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 "AbstractLinAlgPack_MatrixComposite.hpp"
43 #include "AbstractLinAlgPack_SpVectorClass.hpp"
44 #include "AbstractLinAlgPack_VectorStdOps.hpp"
45 #include "AbstractLinAlgPack_VectorMutableBlocked.hpp"
46 #include "AbstractLinAlgPack_AssertOp.hpp"
47 //#include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
48 #include "Teuchos_Workspace.hpp"
49 #include "Teuchos_Assert.hpp"
50 #include "ProfileHackPack_profile_hack.hpp"
51 
52 namespace {
53 
54 // Get an element from a vector (two versions)
55 
56 inline
57 AbstractLinAlgPack::value_type
58 get_element( const AbstractLinAlgPack::Vector& v, AbstractLinAlgPack::index_type i )
59 {
60  return v.get_ele(i);
61 }
62 
63 inline
64 AbstractLinAlgPack::value_type
65 get_element( const AbstractLinAlgPack::SpVectorSlice& v, AbstractLinAlgPack::index_type i )
66 {
68  *ele = v.lookup_element(i);
69  return ele != NULL ? ele->value() : 0.0;
70 }
71 
72 // Get a view of a vector (two versions)
73 
74 inline
76 get_view(
78  ,AbstractLinAlgPack::index_type l
79  ,AbstractLinAlgPack::index_type u
80  )
81 {
82  return v.sub_view(l,u);
83 }
84 
85 inline
87 get_view(
89  ,AbstractLinAlgPack::index_type l
90  ,AbstractLinAlgPack::index_type u
91  )
92 {
93  return Teuchos::rcp( new AbstractLinAlgPack::SpVectorSlice( v(l,u) ) );
94 }
95 
96 // Perform a matrix vector multiplication
97 
98 template<class V>
99 void Vp_StMtV_imp(
100  AbstractLinAlgPack::VectorMutable* y, AbstractLinAlgPack::value_type a
101  ,AbstractLinAlgPack::size_type M_rows, AbstractLinAlgPack::size_type M_cols
104  ,BLAS_Cpp::Transp M_trans
105  ,const V& x, AbstractLinAlgPack::value_type b
106  )
107 {
108  using BLAS_Cpp::no_trans;
109  using BLAS_Cpp::trans;
110  using BLAS_Cpp::rows;
111  using BLAS_Cpp::cols;
112  using AbstractLinAlgPack::dot; // We should not have to do this but some compiles &%!#$
115 
116  AbstractLinAlgPack::Vt_S( y, b ); // Will take care of b == 0.0
117 
118  if( vec_list.size() ) {
119  for( vec_list_t::const_iterator itr = vec_list.begin(); itr != vec_list.end(); ++itr ) {
120  const BLAS_Cpp::Transp
121  op_v_trans = ( M_trans == itr->v_trans_ ? no_trans : trans );
122  const AbstractLinAlgPack::index_type
123  r = ( M_trans == no_trans ? itr->r_l_ : itr->c_l_ ),
124  c = ( M_trans == no_trans ? itr->c_l_ : itr->r_l_ );
125  if( itr->rng_G_.full_range() ) { // op(v)
126  if( op_v_trans == no_trans ) {
127  //
128  // [ y1 ] [ ] [ x1 ]
129  // r:r+n-1 [ y2 ] += a * [ beta * v ] * [ x2 ] c:c
130  // [ y3 ] [ ] [ x3 ]
131  // \______/
132  // c:c
133  // =>
134  //
135  // y(r:r+n-1) += (a * beta * x(c)) * v
136  //
137  AbstractLinAlgPack::Vp_StV( y->sub_view(r,r+itr->v_->dim()-1).get(), a * itr->beta_ * get_element(x,c), *itr->v_ );
138  }
139  else {
140  //
141  // [ y1 ] [ ] [ x1 ]
142  // r:r [ y2 ] += a * [ beta*v' ] * [ x2 ] c:c+n-1
143  // [ y3 ] [ ] [ x3 ]
144  // \_____/
145  // c:c+n-1
146  // =>
147  //
148  // y(r) += a * beta * v'*x(c,c+n-1)
149  //
150 // y->set_ele( r, y->get_ele(r) + a * itr->beta_ * dot( *itr->v_, *get_view(x,c,c+itr->v_->dim()-1) ) );
151  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement the above method in VectorStdOps for Vector,SpVectorSlice!
152  }
153  }
154  else { // op(op(G)*v) or op(v(rng_G))
155  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
156  }
157  }
158  }
159  if( mat_list.size() ) {
160  for( mat_list_t::const_iterator itr = mat_list.begin(); itr != mat_list.end(); ++itr ) {
161  const AbstractLinAlgPack::index_type
162  rl = rows(itr->r_l_,itr->c_l_,M_trans),
163  ru = rows(itr->r_u_,itr->c_u_,M_trans),
164  cl = cols(itr->r_l_,itr->c_l_,M_trans),
165  cu = cols(itr->r_u_,itr->c_u_,M_trans);
166  const BLAS_Cpp::Transp
167  op_P_trans = ( M_trans == itr->P_trans_ ? no_trans : trans ),
168  op_A_trans = ( M_trans == itr->A_trans_ ? no_trans : trans ),
169  op_Q_trans = ( M_trans == itr->Q_trans_ ? no_trans : trans );
170  if( itr->rng_P_.full_range() && itr->rng_Q_.full_range() ) { // op(A)
171  //
172  // [ y1 ] [ ] [ x1 ]
173  // rl:ru [ y2 ] += a * [ alpha * op(op(A)) ] * [ x2 ] cl:cu
174  // [ y3 ] [ ] [ x3 ]
175  // \_______________/
176  // cl:cu
177  // =>
178  //
179  // y(rl:ru) += a * alpha * op(op(A)) * x(cl:cu)
180  //
181  AbstractLinAlgPack::Vp_StMtV( y->sub_view(rl,ru).get(), a * itr->alpha_, *itr->A_, op_A_trans, *get_view(x,cl,cu) );
182  }
183  else {
184  if( itr->A_ == NULL ) { // op(P)
185  TEUCHOS_TEST_FOR_EXCEPT( !( itr->P_.get() && !itr->P_->is_identity() ) );
186  //
187  // [ y1 ] [ ] [ x1 ]
188  // rl:ru [ y2 ] += a * [ alpha * op(op(P)) ] * [ x2 ] cl:cu
189  // [ y3 ] [ ] [ x3 ]
190  // \_______________/
191  // cl:cu
192  // =>
193  //
194  // y(rl:ru) += a * alpha * op(op(P)) * x(cl:cu)
195  //
196 // AbstractLinAlgPack::Vp_StMtV( y->sub_view(rl,ru).get(), a * itr->alpha_, itr->P_, op_P_trans, *get_view(x,cl,cu) );
197  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement the above method properly!
198  }
199  else { // op(P)*op(A)*op(Q) [or some simplification]
200  //
201  // [ y1 ] [ ] [ x1 ]
202  // rl:ru [ y2 ] += a * [ alpha * op(P)*op(A)*op(Q) ] * [ x2 ] cl:cu
203  // [ y3 ] [ ] [ x3 ]
204  // \______________________/
205  // cl:cu
206  // =>
207  //
208  // y(rl:ru) += a * alpha * op(P)*op(A)*op(Q) * x(cl:cu)
209  //
210  if( !itr->rng_P_.full_range() && !itr->rng_Q_.full_range() ) { // op(A)(rng_Q,rng_Q)
212  y->sub_view(rl,ru).get(), a * itr->alpha_
213  ,*itr->A_->sub_view(
214  itr->A_trans_==no_trans ? itr->rng_P_ : itr->rng_Q_
215  ,itr->A_trans_==no_trans ? itr->rng_Q_ : itr->rng_P_ )
216  ,op_A_trans
217  ,*get_view(x,cl,cu)
218  );
219  }
220  else {
221  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
222  }
223  }
224  }
225  }
226  }
227 }
228 
229 // Define a comparison object for comparing SubVectorEntry or SubMatrixEntry
230 // objects for storting by row or by column
231 template<class T>
232 class CompSubEntry {
233 public:
234  enum EByRowCol { BY_ROW, BY_COL };
235  CompSubEntry(enum EByRowCol by_row_col)
236  : by_row_col_(by_row_col)
237  {}
238  bool operator()( const T& e1, const T& e2 )
239  {
240  return
241  ( by_row_col_ == BY_ROW
242  ? e1.r_l_ < e2.r_l_
243  : e1.c_l_ < e2.c_l_
244  );
245  }
246 private:
247  EByRowCol by_row_col_;
248  CompSubEntry(); // not defined and not to be called
249 }; // end class CompSubVectorEntry
250 
251 
252 } // end namespace
253 
254 namespace AbstractLinAlgPack {
255 
257 {
258  reinitialize(rows,cols);
259 }
260 
261 void MatrixComposite::reinitialize( size_type rows, size_type cols )
262 {
263  namespace rcp = MemMngPack;
264 
265  fully_constructed_ = true;
266  rows_ = rows;
267  cols_ = cols;
268  if(matrix_list_.size()) {
269  matrix_list_.erase(matrix_list_.begin(),matrix_list_.end());
270  }
271  if(vector_list_.size()) {
272  vector_list_.erase(vector_list_.begin(),vector_list_.end());
273  }
274  space_rows_ = Teuchos::null;
275  space_cols_ = Teuchos::null;
276 }
277 
279  size_type row_offset
280  ,size_type col_offset
281  ,value_type beta
282  ,const GenPermMatrixSlice *G
283  ,const release_resource_ptr_t &G_release
284  ,BLAS_Cpp::Transp G_trans
285  ,const Vector *v
286  ,const release_resource_ptr_t &v_release
287  ,BLAS_Cpp::Transp v_trans
288  )
289 {
290  fully_constructed_ = false;
291  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
292 }
293 
295  size_type row_offset
296  ,size_type col_offset
297  ,value_type beta
298  ,const Range1D &rng_G
299  ,const Vector *v
300  ,const release_resource_ptr_t &v_release
301  ,BLAS_Cpp::Transp v_trans
302  )
303 {
304  fully_constructed_ = false;
305  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
306 }
307 
309  size_type row_offset
310  ,size_type col_offset
311  ,value_type beta
312  ,const Vector *v
313  ,const release_resource_ptr_t &v_release
314  ,BLAS_Cpp::Transp v_trans
315  )
316 {
317  namespace rcp = MemMngPack;
318 
319  TEUCHOS_TEST_FOR_EXCEPT( !( beta != 0.0 ) );
320  TEUCHOS_TEST_FOR_EXCEPT( !( v != NULL ) );
321  fully_constructed_ = false;
322  if( v_trans == BLAS_Cpp::no_trans ) {
323  TEUCHOS_TEST_FOR_EXCEPT( !( row_offset + v->dim() <= rows_ ) );
324  TEUCHOS_TEST_FOR_EXCEPT( !( col_offset + 1 <= cols_ ) );
325  }
326  else {
327  TEUCHOS_TEST_FOR_EXCEPT( !( row_offset + 1 <= rows_ ) );
328  TEUCHOS_TEST_FOR_EXCEPT( !( col_offset + v->dim() <= cols_ ) );
329  }
330 
331  vector_list_.push_back(
333  row_offset+1,col_offset+1,beta,Range1D()
334  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
335  ,v,v_release,v_trans ) );
336 }
337 
338 void MatrixComposite::remove_vector( vector_list_t::iterator itr )
339 {
340  fully_constructed_ = false;
341  vector_list_.erase(itr);
342 }
343 
345  size_type row_offset
346  ,size_type col_offset
347  ,value_type alpha
348  ,const GenPermMatrixSlice *P
349  ,const release_resource_ptr_t &P_release
350  ,BLAS_Cpp::Transp P_trans
351  ,const MatrixOp *A
352  ,const release_resource_ptr_t &A_release
353  ,BLAS_Cpp::Transp A_trans
354  ,const GenPermMatrixSlice *Q
355  ,const release_resource_ptr_t &Q_release
356  ,BLAS_Cpp::Transp Q_trans
357  )
358 {
359  fully_constructed_ = false;
360  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
361 }
362 
364  size_type row_offset
365  ,size_type col_offset
366  ,value_type alpha
367  ,const Range1D &rng_P
368  ,const MatrixOp *A
369  ,const release_resource_ptr_t &A_release
370  ,BLAS_Cpp::Transp A_trans
371  ,const GenPermMatrixSlice *Q
372  ,const release_resource_ptr_t &Q_release
373  ,BLAS_Cpp::Transp Q_trans
374  )
375 {
376  fully_constructed_ = false;
377  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
378 }
379 
381  size_type row_offset
382  ,size_type col_offset
383  ,value_type alpha
384  ,const GenPermMatrixSlice *P
385  ,const release_resource_ptr_t &P_release
386  ,BLAS_Cpp::Transp P_trans
387  ,const MatrixOp *A
388  ,const release_resource_ptr_t &A_release
389  ,BLAS_Cpp::Transp A_trans
390  ,const Range1D &rng_Q
391  )
392 {
393  fully_constructed_ = false;
394  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish!
395 }
396 
398  size_type row_offset
399  ,size_type col_offset
400  ,value_type alpha
401  ,const Range1D &rng_P_in
402  ,const MatrixOp *A
403  ,const release_resource_ptr_t &A_release
404  ,BLAS_Cpp::Transp A_trans
405  ,const Range1D &rng_Q_in
406  )
407 {
408  namespace rcp = MemMngPack;
409  using BLAS_Cpp::rows;
410  using BLAS_Cpp::cols;
411  using RangePack::full_range;
412 
414  alpha == 0.0, std::invalid_argument
415  ,"MatrixComposite::add_matrix(...) : Error!" );
417  A == NULL, std::invalid_argument
418  ,"MatrixComposite::add_matrix(...) : Error!" );
419 
420  const size_type
421  A_rows = A->rows(),
422  A_cols = A->cols(),
423  opA_rows = rows(A_rows,A_cols,A_trans),
424  opA_cols = cols(A_rows,A_cols,A_trans);
425  const Range1D
426  rng_P = full_range(rng_P_in,1,opA_rows),
427  rng_Q = full_range(rng_Q_in,1,opA_cols);
428  const size_type
429  opPopAopQ_rows = rng_P.size(),
430  opPopAopQ_cols = rng_Q.size();
431 
433  row_offset + opPopAopQ_rows > rows_, std::invalid_argument
434  ,"MatrixComposite::add_matrix(...) : Error!" );
436  col_offset + opPopAopQ_cols > cols_, std::invalid_argument
437  ,"MatrixComposite::add_matrix(...) : Error!" );
438 
439  fully_constructed_ = false;
440 
441  matrix_list_.push_back(
443  row_offset+1,row_offset+opPopAopQ_rows
444  ,col_offset+1,col_offset+opPopAopQ_cols
445  ,alpha
446  ,rng_P
447  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
448  ,A,A_release,A_trans
449  ,rng_Q
450  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
451  )
452  );
453 
454  // ToDo: We could create a sub-view of the matrix using rng_P and rng_Q
455  // and then store the sub-view in the data structure? However, this
456  // would confuse an external client who wanted to iterate through
457  // the matrix list using the public iterators.
458 }
459 
461  size_type row_offset
462  ,size_type col_offset
463  ,value_type alpha
464  ,const MatrixOp *A
465  ,const release_resource_ptr_t &A_release
466  ,BLAS_Cpp::Transp A_trans
467  )
468 {
469  namespace rcp = MemMngPack;
470  using BLAS_Cpp::rows;
471  using BLAS_Cpp::cols;
472 
474  alpha == 0.0, std::invalid_argument
475  ,"MatrixComposite::add_matrix(...) : Error!" );
477  A == NULL, std::invalid_argument
478  ,"MatrixComposite::add_matrix(...) : Error!" );
479 
480  const size_type
481  A_rows = A->rows(),
482  A_cols = A->cols(),
483  opA_rows = rows(A_rows,A_cols,A_trans),
484  opA_cols = cols(A_rows,A_cols,A_trans);
485 
487  row_offset + opA_rows > rows_, std::invalid_argument
488  ,"MatrixComposite::add_matrix(...) : Error!" );
490  col_offset + opA_cols > cols_, std::invalid_argument
491  ,"MatrixComposite::add_matrix(...) : Error!" );
492 
493  fully_constructed_ = false;
494 
495  matrix_list_.push_back(
497  row_offset+1,row_offset+opA_rows
498  ,col_offset+1,col_offset+opA_cols
499  ,alpha
500  ,Range1D()
501  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
502  ,A,A_release,A_trans
503  ,Range1D()
504  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
505  )
506  );
507 }
508 
510  size_type row_offset
511  ,size_type col_offset
512  ,value_type alpha
513  ,const GenPermMatrixSlice *P
514  ,const release_resource_ptr_t &P_release
515  ,BLAS_Cpp::Transp P_trans
516  )
517 {
518  namespace rcp = MemMngPack;
519  using BLAS_Cpp::rows;
520  using BLAS_Cpp::cols;
521 
522  TEUCHOS_TEST_FOR_EXCEPT( !( alpha != 0.0 ) );
523  TEUCHOS_TEST_FOR_EXCEPT( !( P != NULL ) );
524 
525  fully_constructed_ = false;
526 
527  const size_type
528  P_rows = P->rows(),
529  P_cols = P->cols(),
530  opP_rows = rows(P_rows,P_cols,P_trans),
531  opP_cols = cols(P_rows,P_cols,P_trans);
532 
533  TEUCHOS_TEST_FOR_EXCEPT( !( row_offset + opP_rows <= rows_ ) );
534  TEUCHOS_TEST_FOR_EXCEPT( !( col_offset + opP_cols <= cols_ ) );
535 
536  matrix_list_.push_back(
538  row_offset+1,row_offset+opP_rows,col_offset+1,col_offset+opP_cols,alpha
540  ,Teuchos::rcp(new GenPermMatrixSlice(*P)),P_release,P_trans
541  ,NULL,Teuchos::null,BLAS_Cpp::no_trans
542  ,Range1D()
543  ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
544  )
545  );
546 }
547 
548 void MatrixComposite::remove_matrix( matrix_list_t::iterator itr )
549 {
550  fully_constructed_ = false;
551  matrix_list_.erase(itr);
552 }
553 
555  const VectorSpace::space_ptr_t& space_cols
556  ,const VectorSpace::space_ptr_t& space_rows
557  )
558 {
560  !space_cols.get(), std::invalid_argument
561  ,"MatrixComposite::finish_construction(...): Error, space_cols.get() can not be NULL" );
563  !space_rows.get(), std::invalid_argument
564  ,"MatrixComposite::finish_construction(...): Error, space_rows.get() can not be NULL" );
566  space_cols->dim() != rows_, std::invalid_argument
567  ,"MatrixComposite::finish_construction(...): Error, space_colss->dim() = " << space_cols->dim()
568  << " != rows = " << rows_ << " where cols was passed to this->reinitialize(...)" );
570  space_rows->dim() != cols_, std::invalid_argument
571  ,"MatrixComposite::finish_construction(...): Error, space_rows->dim() = " << space_rows->dim()
572  << " != cols = " << cols_ << " where cols was passed to this->reinitialize(...)" );
573 
574  space_rows_ = space_rows;
575  space_cols_ = space_cols;
576  fully_constructed_ = true;
577 }
578 
579 // Member access
580 
582 {
583  return vector_list_.size();
584 }
585 
586 MatrixComposite::vector_list_t::iterator
588 {
589  return vector_list_.begin();
590 }
591 
592 MatrixComposite::vector_list_t::iterator
594 {
595  return vector_list_.end();
596 }
597 
598 MatrixComposite::vector_list_t::const_iterator
600 {
601  return vector_list_.begin();
602 }
603 
604 MatrixComposite::vector_list_t::const_iterator
606 {
607  return vector_list_.end();
608 }
609 
611 {
612  return matrix_list_.size();
613 }
614 
615 MatrixComposite::matrix_list_t::iterator
617 {
618  return matrix_list_.begin();
619 }
620 
621 MatrixComposite::matrix_list_t::iterator
623 {
624  return matrix_list_.end();
625 }
626 
627 MatrixComposite::matrix_list_t::const_iterator
629 {
630  return matrix_list_.begin();
631 }
632 
633 MatrixComposite::matrix_list_t::const_iterator
635 {
636  return matrix_list_.end();
637 }
638 
639 // Overridden from Matrix
640 
641 size_type MatrixComposite::rows() const
642 {
643  return fully_constructed_ ? rows_ : 0;
644 }
645 
646 size_type MatrixComposite::cols() const
647 {
648  return fully_constructed_ ? cols_ : 0;
649 }
650 
651 size_type MatrixComposite::nz() const
652 {
653  if(fully_constructed_) {
654  size_type nz = 0;
655  {for(matrix_list_t::const_iterator itr = matrix_list_.begin(); itr != matrix_list_.end(); ++itr) {
656  if( itr->A_ != NULL )
657  nz += itr->A_->nz();
658  else
659  nz += (*itr->P_).nz();
660  }}
661  {for(vector_list_t::const_iterator itr = vector_list_.begin(); itr != vector_list_.end(); ++itr) {
662  nz += itr->v_->nz();
663  }}
664  // ToDo: Update the above code to consider DMatrixSlice and Range1D views
665  return nz;
666  }
667  return 0;
668 }
669 
670 // Overridden from MatrixOp
671 
673 {
674  assert_fully_constructed();
675  return *space_rows_;
676 }
677 
678 const VectorSpace& MatrixComposite::space_cols() const
679 {
680  assert_fully_constructed();
681  return *space_cols_;
682 }
683 
684 MatrixOp::mat_ptr_t
685 MatrixComposite::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
686 {
687  assert_fully_constructed();
688  // For this initial implementation we will just look through the list of matrices
689  // and find an exact match. If we don't find it we will just return the result
690  // from the default implementation of this method.
691 
692  // ToDo: Implement!
693 
694  // Just return the default sub-view
695  return MatrixOp::sub_view(row_rng,col_rng);
696 }
697 
699  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
700  , const Vector& x, value_type b
701  ) const
702 {
703 #ifdef PROFILE_HACK_ENABLED
704  ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StMtV(...DVectorSlice...)" );
705 #endif
706  assert_fully_constructed();
707  AbstractLinAlgPack::Vp_MtV_assert_compatibility( y, *this, M_trans, x );
708  Vp_StMtV_imp(y,a,rows_,cols_,matrix_list_,vector_list_,M_trans,x,b);
709 }
710 
712  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
713  , const SpVectorSlice& x, value_type b
714  ) const
715 {
716 #ifdef PROFILE_HACK_ENABLED
717  ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StMtV(...SpVectorSlice...)" );
718 #endif
719  assert_fully_constructed();
720  AbstractLinAlgPack::Vp_MtV_assert_compatibility( y, *this, M_trans, x );
721  Vp_StMtV_imp(y,a,rows_,cols_,matrix_list_,vector_list_,M_trans,x,b);
722 }
723 
725  VectorMutable* y, value_type a
726  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
727  , BLAS_Cpp::Transp M_trans
728  , const Vector& x, value_type b
729  ) const
730 {
731 #ifdef PROFILE_HACK_ENABLED
732  ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StPtMtV(...DVectorSlice...)" );
733 #endif
734  assert_fully_constructed();
735  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement when needed!
736 }
737 
739  VectorMutable* y, value_type a
740  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
741  , BLAS_Cpp::Transp M_trans
742  , const SpVectorSlice& x, value_type b
743  ) const
744 {
745 #ifdef PROFILE_HACK_ENABLED
746  ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StPtMtV(...SpVectorSlice...)" );
747 #endif
748  assert_fully_constructed();
749  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement when needed!
750 }
751 
752 // private
753 
754 void MatrixComposite::assert_fully_constructed() const
755 {
756  const bool fully_constructed = fully_constructed_;
758  !fully_constructed, std::logic_error
759  ,"MatrixComposite::assert_fully_constructed() : Error, not fully constructed!");
760 }
761 
762 } // end namespace AbstractLinAlgPack
void remove_vector(vector_list_t::iterator itr)
Remove a sub-vector.
virtual vec_ptr_t sub_view(const Range1D &rng) const
Create an abstract view of a vector object .
virtual vec_mut_ptr_t sub_view(const Range1D &rng)
Create a mutable abstract view of a vector object.
std::deque< SubMatrixEntry > matrix_list_t
Warning! This could be changed to some other STL container!
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
virtual mat_ptr_t sub_view(const Range1D &row_rng, const Range1D &col_rng) const
Create a transient constant sub-matrix view of this matrix (if supported).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index size() const
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
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).
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
mat_ptr_t sub_view(const Range1D &row_rng, const Range1D &col_rng) const
T * get() const
virtual void finish_construction(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows)
Call to finish the construction process.
std::deque< SubVectorEntry > vector_list_t
Warning! This could be changed to some other STL container!
virtual size_type cols() const
Return the number of columns in the matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract interface for objects that represent a space for mutable coordinate vectors.
void Vp_StPtMtV(VectorMutable *vs_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
const VectorSpace & space_cols() const
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
Base class for all matrices that support basic matrix operations.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
MatrixComposite(size_type rows=0, size_type cols=0)
Construct.
static const Range1D Invalid
void reinitialize(size_type rows=0, size_type cols=0)
Initialize a sized (on unsized) zero matrix to start with.
virtual index_type dim() const
Return the dimension of this vector.
virtual value_type get_ele(index_type i) const
Fetch an element in the vector.
Abstract interface for mutable coordinate vectors {abstract}.
void Vp_MtV_assert_compatibility(VectorMutable *v_lhs, const MatrixOp &m_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs += op(m_rhs1) * v_rhs2
virtual size_type rows() const
Return the number of rows in the matrix.
void remove_matrix(matrix_list_t::iterator itr)
Remove a sub-matrix.
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void add_vector(size_type row_offset, size_type col_offset, value_type beta, const GenPermMatrixSlice *G, const release_resource_ptr_t &G_release, BLAS_Cpp::Transp G_trans, const Vector *v, const release_resource_ptr_t &v_release, BLAS_Cpp::Transp v_trans)
Add a sub-vector beta*op(op(G)*v).
const VectorSpace & space_rows() const
Concrete matrix type to represent general permutation (mapping) matrices.
friend 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)
void Vp_StMtV(VectorMutable *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const