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_MultiVectorMutableCols.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_MultiVectorMutableCols.hpp"
43 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
44 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
45 #include "AbstractLinAlgPack_SpVectorClass.hpp"
46 #include "Teuchos_dyn_cast.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 namespace AbstractLinAlgPack {
50 
51 // Constructors/Initializers
52 
54 {}
55 
57  const Teuchos::RCP<const VectorSpace> &space_cols
58  ,const Teuchos::RCP<const VectorSpace> &space_rows
59  ,Teuchos::RCP<VectorMutable> col_vecs[]
60  )
61 {
62  this->initialize(space_cols,space_rows,col_vecs);
63 }
64 
66  const Teuchos::RCP<const VectorSpace> &space_cols
67  ,const Teuchos::RCP<const VectorSpace> &space_rows
68  ,Teuchos::RCP<VectorMutable> col_vecs[]
69  )
70 {
71  space_cols_ = space_cols;
72  space_rows_ = space_rows;
73  const size_type num_cols = space_rows->dim();
74  col_vecs_.resize(num_cols);
75  if(col_vecs) {
76  for( size_type j = 1; j <= num_cols; ++j )
77  col_vecs_[j-1] = col_vecs[j-1];
78  }
79  else {
80  for( size_type j = 1; j <= num_cols; ++j )
81  col_vecs_[j-1] = space_cols->create_member();
82  }
83 }
84 
86 {
87  col_vecs_.resize(0);
88  space_cols_ = Teuchos::null;
89  space_rows_ = Teuchos::null;
90 }
91 
92 // Overridden from MatrixBase
93 
95 {
96  return space_cols_.get() ? space_cols_->dim() : 0;
97 }
98 
100 {
101  return space_rows_.get() ? space_rows_->dim() : 0;
102 }
103 
104 // Overridden from MatrixOp
105 
107 {
108  return *space_cols_;
109 }
110 
111 const VectorSpace& MultiVectorMutableCols::space_rows() const
112 {
113  return *space_rows_;
114 }
115 
117 {
118  for( size_type j = 1; j <= col_vecs_.size(); ++j )
119  col_vecs_[j-1]->zero();
120 }
121 
122 void MultiVectorMutableCols::Mt_S( value_type alpha )
123 {
124  for( size_type j = 1; j <= col_vecs_.size(); ++j )
125  LinAlgOpPack::Vt_S(col_vecs_[j-1].get(),alpha);
126 }
127 
128 MatrixOp&
130 {
131  const MultiVector
132  *mv_rhs = dynamic_cast<const MultiVector*>(&mwo_rhs);
133  if(!mv_rhs)
134  return MatrixOp::operator=(mwo_rhs);
135  for( size_type j = 1; j <= col_vecs_.size(); ++j )
136  *col_vecs_[j-1] = *mv_rhs->col(j);
137  return *this;
138 }
139 
140 MatrixOp::mat_mut_ptr_t
142 {
143  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
144  return Teuchos::null;
145 }
146 
147 MatrixOp::mat_ptr_t
149 {
150  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
151  return Teuchos::null;
152 }
153 
155  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
156  ,const Vector& x, value_type b
157  ) const
158 {
160 
161  // y = b*y
162  LinAlgOpPack::Vt_S(y,b);
163 
164  if( M_trans == BLAS_Cpp::no_trans ) {
165  //
166  // y += a*M*x
167  //
168  // =>
169  //
170  // y += a * x(1) * M(:,1) + a * x(2) * M(:,2) + ...
171  //
172  for( size_type j = 1; j <= col_vecs_.size(); ++j )
173  LinAlgOpPack::Vp_StV( y, a * x.get_ele(j), *col_vecs_[j-1] );
174  }
175  else {
176  //
177  // y += a*M'*x
178  //
179  // =>
180  //
181  // y(1) += a M(:,1)*x
182  // y(2) += a M(:,2)*x
183  // ...
184  //
185  for( size_type j = 1; j <= col_vecs_.size(); ++j )
186  y->set_ele(
187  j
188  ,y->get_ele(j) + a * dot(*col_vecs_[j-1],x)
189  );
190  }
191 }
192 
194  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
195  ,const SpVectorSlice& x, value_type b
196  ) const
197 {
199 
200  // y = b*y
201  LinAlgOpPack::Vt_S(y,b);
202 
203  if( M_trans == BLAS_Cpp::no_trans ) {
204  //
205  // y += a*M*x
206  //
207  // =>
208  //
209  // y += a * x(1) * M(:,1) + a * x(2) * M(:,2) + ...
210  //
212  for( SpVectorSlice::const_iterator itr = x.begin(); itr != x.end(); ++itr ) {
213  const size_type j = itr->index() + o;
214  LinAlgOpPack::Vp_StV( y, a * itr->value(), *col_vecs_[j-1] );
215  }
216  }
217  else {
218  //
219  // y += a*M'*x
220  //
221  // =>
222  //
223  // y(1) += a M(:,1)*x
224  // y(2) += a M(:,2)*x
225  // ...
226  //
227  for( size_type j = 1; j <= col_vecs_.size(); ++j )
228  y->set_ele(
229  j
230  ,y->get_ele(j) + a * dot(*col_vecs_[j-1],x)
231  );
232  }
233 }
234 
236  BLAS_Cpp::Transp M_trans, value_type alpha
237  , value_type beta, MatrixSymOp* sym_lhs ) const
238 {
239  using LinAlgOpPack::dot;
241  *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(sym_lhs);
242  if(!symwo_gms_lhs) {
243  return MatrixOp::syrk(M_trans,alpha,beta,sym_lhs); // Boot it
244  }
246  const int num_vecs = this->col_vecs_.size();
248  num_vecs != DMatrixSliceSym().rows(), std::logic_error
249  ,"MultiVectorMutableCols::syrk(...) : Error, sizes do not match up!" );
250  // Fill the upper or lower triangular region.
251  if( DMatrixSliceSym().uplo() == BLAS_Cpp::upper ) {
252  for( int i = 1; i <= num_vecs; ++i ) {
253  for( int j = i; j <= num_vecs; ++j ) { // Upper triangle!
254  DMatrixSliceSym().gms()(i,j) = beta * DMatrixSliceSym().gms()(i,j) + alpha * dot(*col_vecs_[i-1],*col_vecs_[j-1]);
255  }
256  }
257  }
258  else {
259  for( int i = 1; i <= num_vecs; ++i ) {
260  for( int j = 1; j <= i; ++j ) { // Lower triangle!
261  DMatrixSliceSym().gms()(i,j) = beta * DMatrixSliceSym().gms()(i,j) + alpha * dot(*col_vecs_[i-1],*col_vecs_[j-1]);
262  }
263  }
264  }
265  return true;
266 }
267 
268 // Overridden from MultiVector
269 
272 {
273  return COL_ACCESS;
274 }
275 
276 // Overridden from MultiVectorMutable
277 
280 {
281  TEUCHOS_TEST_FOR_EXCEPTION( !( 1 <= j && j <= col_vecs_.size() ), std::logic_error, "Error!" );
282  return col_vecs_[j-1];
283 }
284 
287 {
288  return Teuchos::null;
289 }
290 
293 {
294  return Teuchos::null;
295 }
296 
298 MultiVectorMutableCols::mv_sub_view(const Range1D& row_rng, const Range1D& col_rng)
299 {
300 #ifdef TEUCHOS_DEBUG
301  const size_type rows = this->rows();
303  !( row_rng.full_range() || (row_rng.lbound() == 1 && row_rng.ubound() == rows) )
304  ,std::logic_error, "Error, can't handle subrange on vectors yet!" );
305 #endif
306  return Teuchos::rcp(
308  space_cols_,space_rows_->sub_space(col_rng),&col_vecs_[col_rng.lbound()-1]
309  ) );
310 }
311 
312 } // end namespace AbstractLinAlgPack
difference_type offset() const
Return the offset for the indexes (ith real index = begin()[i-1]->index() + offset()# ...
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
bool syrk(BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
Index ubound() const
virtual MatrixOp & operator=(const MatrixOp &mwo_rhs)
M_lhs = mwo_rhs : Virtual assignment operator.
Abstract interface that allows the extraction of a non-const DenseLinAlgPack::DMatrixSliceSym view of...
virtual void set_ele(index_type i, value_type val)
Set a specific element of a vector.
Interface adding operations specific for a symmetric matrix {abstract}.
vec_mut_ptr_t diag(int k)
Returns return.get() == NULL
void initialize(const Teuchos::RCP< const VectorSpace > &space_cols, const Teuchos::RCP< const VectorSpace > &space_rows, Teuchos::RCP< VectorMutable > col_vecs[]=NULL)
Initialize given the spaces for the columns and rows and possibly the column vectors.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const VectorSpace & space_cols() const
Abstract interface for objects that represent a space for mutable coordinate vectors.
Helper class type that simplifies the usage of the MatrixSymOpGetGMSSymMutable interface for clients...
const VectorSpace & space_rows() const
bool full_range() const
access_by_t access_by() const
Returns return & COL_ACCESS == true only.
Base class for all matrices that support basic matrix operations.
Interface for a collection of non-mutable vectors (multi-vector, matrix).
multi_vec_mut_ptr_t mv_sub_view(const Range1D &row_rng, const Range1D &col_rng)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Index lbound() const
virtual value_type get_ele(index_type i) const
Fetch an element in the vector.
Abstract interface for mutable coordinate vectors {abstract}.
virtual vec_ptr_t col(index_type j) const =0
Get a non-mutable column vector.
vec_mut_ptr_t row(index_type i)
Returns return.get() == NULL
friend void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Transp
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)