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_MatrixSymDiagonalStd.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymDiagStd.hpp"
45 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
46 #include "DenseLinAlgPack_DVectorOp.hpp"
47 #include "DenseLinAlgPack_DMatrixOp.hpp"
48 #include "DenseLinAlgPack_AssertOp.hpp"
49 
50 namespace AbstractLinAlgPack {
51 
53 {}
54 
55 DVectorSlice MatrixSymDiagStd::diag()
56 {
57  return diag_();
58 }
59 
60 const DVectorSlice MatrixSymDiagStd::diag() const
61 {
62  return diag_();
63 }
64 
65 // Overridden from MatrixSymInitDiag
66 
67 void MatrixSymDiagStd::init_identity( size_type n, value_type alpha = 1.0 )
68 {
69  diag_.resize(n);
70  if(n)
71  diag_ = alpha;
72 }
73 
74 void MatrixSymDiagStd::init_diagonal( const DVectorSlice& diag )
75 {
76  diag_ = diag;
77 }
78 
79 // Overridden from Matrix
80 
82 {
83  return diag_.size();
84 }
85 
86 // Overridden from MatrixOp
87 
89  DMatrixSlice* gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
90 {
92  // Assert that the dimensions match up.
93  DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
94  , rows(), cols(), trans_rhs );
95  // Add to the diagonal
96  Vp_StV( &gms_lhs->diag(), alpha, diag_ );
97 }
98 
100  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
101  , const DVectorSlice& vs_rhs2, value_type beta) const
102 {
103  const DVectorSlice diag = this->diag();
104  size_type n = diag.size();
105 
106  //
107  // y = b*y + a * op(A) * x
108  //
109  DenseLinAlgPack::Vp_MtV_assert_sizes(
110  vs_lhs->size(), n, n, trans_rhs1, vs_rhs2.size() );
111  //
112  // A is symmetric and diagonal A = diag(diag) so:
113  //
114  // y(j) += a * diag(j) * x(j), for j = 1...n
115  //
116  if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
117  // Optimized implementation
118  const value_type
119  *d_itr = diag.raw_ptr(),
120  *x_itr = vs_rhs2.raw_ptr();
121  value_type
122  *y_itr = vs_lhs->raw_ptr(),
123  *y_end = y_itr + vs_lhs->size();
124 
125  if( beta == 0.0 ) {
126  while( y_itr != y_end )
127  *y_itr++ = alpha * (*d_itr++) * (*x_itr++);
128  }
129  else if( beta == 1.0 ) {
130  while( y_itr != y_end )
131  *y_itr++ += alpha * (*d_itr++) * (*x_itr++);
132  }
133  else {
134  for( ; y_itr != y_end; ++y_itr )
135  *y_itr = beta * (*y_itr) + alpha * (*d_itr++) * (*x_itr++);
136  }
137  }
138  else {
139  // Generic implementation
140  DVectorSlice::const_iterator
141  d_itr = diag.begin(),
142  x_itr = vs_rhs2.begin();
143  DVectorSlice::iterator
144  y_itr = vs_lhs->begin(),
145  y_end = vs_lhs->end();
146  for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
147 #ifdef LINALGPACK_CHECK_RANGE
148  TEUCHOS_TEST_FOR_EXCEPT( !( d_itr < diag.end() ) );
149  TEUCHOS_TEST_FOR_EXCEPT( !( x_itr < vs_rhs2.end() ) );
150  TEUCHOS_TEST_FOR_EXCEPT( !( y_itr < vs_lhs->end() ) );
151 #endif
152  *y_itr = beta * (*y_itr) + alpha * (*d_itr) * (*x_itr);
153  }
154  }
155 }
156 
158  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
159  , const SpVectorSlice& sv_rhs2, value_type beta) const
160 {
161  const DVectorSlice diag = this->diag();
162  size_type n = diag.size();
163 
164  DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
165  , n, n, trans_rhs1, sv_rhs2.size() );
166  //
167  // y = b*y + a * op(A) * x
168  //
169  DenseLinAlgPack::Vt_S(vs_lhs,beta); // y = b * y
170  //
171  // A is symmetric and diagonal A = diag(diag) so:
172  //
173  // y(j) += a * diag(j) * x(j), for j = 1...n
174  //
175  // x is sparse so take account of this.
176 
177  for( SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
178  ; x_itr != sv_rhs2.end()
179  ; ++x_itr )
180  {
181  (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
182  += alpha * diag(x_itr->indice() + sv_rhs2.offset()) * x_itr->value();
183  // Note: The indice x(i) invocations are ranged check
184  // if this is compiled into the code.
185  }
186 }
187 
188 // Overridden from MatrixWithOpFactorized
189 
191  DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
192  , const DVectorSlice& vs_rhs2) const
193 {
194  const DVectorSlice diag = this->diag();
195  size_type n = diag.size();
196 
197  // y = inv(op(A)) * x
198  //
199  // A is symmetric and diagonal (A = diag(diag)) so:
200  //
201  // y(j) = x(j) / diag(j), for j = 1...n
202 
203  DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
204  , n, n, trans_rhs1, vs_rhs2.size() );
205 
206  if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
207  // Optimized implementation
208  const value_type
209  *d_itr = diag.raw_ptr(),
210  *x_itr = vs_rhs2.raw_ptr();
211  value_type
212  *y_itr = vs_lhs->raw_ptr(),
213  *y_end = y_itr + vs_lhs->size();
214  while( y_itr != y_end )
215  *y_itr++ = (*x_itr++) / (*d_itr++);
216  }
217  else {
218  // Generic implementation
219  DVectorSlice::const_iterator
220  d_itr = diag.begin(),
221  x_itr = vs_rhs2.begin();
222  DVectorSlice::iterator
223  y_itr = vs_lhs->begin(),
224  y_end = vs_lhs->end();
225  for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
226  TEUCHOS_TEST_FOR_EXCEPT( !( d_itr < diag.end() ) );
227  TEUCHOS_TEST_FOR_EXCEPT( !( x_itr < vs_rhs2.end() ) );
228  TEUCHOS_TEST_FOR_EXCEPT( !( y_itr < vs_lhs->end() ) );
229  *y_itr = (*x_itr)/(*d_itr);
230  }
231  }
232 }
233 
235  DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
236  , const SpVectorSlice& sv_rhs2) const
237 {
238  const DVectorSlice diag = this->diag();
239  size_type n = diag.size();
240 
241  // y = inv(op(A)) * x
242  //
243  // A is symmetric and diagonal A = diag(diag) so:
244  //
245  // y(j) = x(j) / diag(j), for j = 1...n
246  //
247  // x is sparse so take account of this.
248 
249  DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
250  , n, n, trans_rhs1, sv_rhs2.size() );
251 
252  for( SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
253  ; x_itr != sv_rhs2.end()
254  ; ++x_itr )
255  {
256  (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
257  = x_itr->value() / diag(x_itr->indice() + sv_rhs2.offset());
258  // Note: The indice x(i) invocations are ranged check
259  // if this is compiled into the code.
260  }
261 }
262 
263 } // end namespace AbstractLinAlgPack
264 
265 #endif // 0
size_type cols() const
Returns this->rows()
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
MatrixSymDiagStd(const VectorSpace::vec_mut_ptr_t &diag=Teuchos::null, bool unique=true)
Calls this->initialize().
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
bool Mp_StM(MatrixOp *g_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
Add to a mutable matrix lhs.
size_t size_type
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void init_identity(const VectorSpace &space_diag, value_type alpha)
size_type rows() const
Returns 0 if not initalized (this->diag() == NULL).
Transp
VectorMutable & diag()
Give non-const access to the diagonal vector.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)