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_MatrixSymDiagSparse.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 <assert.h>
43 
44 #include <fstream> // For debugging only
45 #include <limits>
46 
47 #include "AbstractLinAlgPack_MatrixSymDiagSparse.hpp"
48 #include "AbstractLinAlgPack_SpVectorClass.hpp"
49 #include "AbstractLinAlgPack_EtaVector.hpp"
50 #include "AbstractLinAlgPack_AssertOp.hpp"
51 #include "AbstractLinAlgPack_SpVectorOut.hpp"
52 #include "DenseLinAlgPack_DVectorClass.hpp"
53 #include "DenseLinAlgPack_DMatrixAssign.hpp"
54 #include "DenseLinAlgPack_DMatrixClass.hpp"
55 #include "DenseLinAlgPack_DMatrixOp.hpp"
56 #include "DenseLinAlgPack_assert_print_nan_inf.hpp"
57 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
58 #include "Teuchos_Assert.hpp"
59 
60 namespace {
61 template< class T >
62 inline
63 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
64 } // end namespace
65 
66 namespace AbstractLinAlgPack {
67 
69  : num_updates_at_once_(0) // Flag that it is to be determined internally.
70 {}
71 
72 // Overridden from MatrixBase
73 
74 size_type MatrixSymDiagSparse::rows() const
75 {
76  return diag().dim();
77 }
78 
79 // Overridden from MatrixOp
80 
81 std::ostream& MatrixSymDiagSparse::output(std::ostream& out) const
82 {
83  out << "*** Sparse diagonal matrix ***\n"
84  << "diag =\n" << diag();
85  return out;
86 }
87 
88 // Overridden from MatrixOpSerial
89 
90 void MatrixSymDiagSparse::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha
91  , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta) const
92 {
93  const SpVectorSlice &diag = this->diag();
94 
95  size_type n = diag.dim();
96 
97  // Assert that the dimensions of the aruments match up and if not
98  // then throw an excption.
99  DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->dim(), n, n, trans_rhs1, vs_rhs2.dim() );
100 
101  // y = b*y + a * op(A) * x
102  //
103  // A is symmetric and diagonal A = diag(diag) so:
104  //
105  // y(j) = b*y(j) + a * diag(j) * x(j), for j = 1...n
106 
107  for( SpVectorSlice::const_iterator d_itr = diag.begin(); d_itr != diag.end(); ++d_itr )
108  {
109  const size_type i = d_itr->index();
110  (*vs_lhs)(i) = beta * (*vs_lhs)(i) + alpha * d_itr->value() * vs_rhs2(i);
111  }
112 }
113 
114 // Overridden from MatrixSymOpSerial
115 
117  DMatrixSliceSym* B, value_type alpha
118  ,EMatRhsPlaceHolder dummy_place_holder
119  ,const MatrixOpSerial& A, BLAS_Cpp::Transp A_trans
120  ,value_type b
121  ) const
122 {
123  using BLAS_Cpp::rows;
124  using BLAS_Cpp::cols;
125  using BLAS_Cpp::trans_not;
126 
127  using DenseLinAlgPack::nonconst_tri_ele;
129  using DenseLinAlgPack::syrk;
130  using DenseLinAlgPack::assert_print_nan_inf;
131 
132  using LinAlgOpPack::V_MtV;
133 
134  typedef EtaVector eta;
135 
136  // Assert the size matching of M * op(A)
137  DenseLinAlgPack::MtV_assert_sizes(
138  this->rows(), this->cols(), BLAS_Cpp::no_trans
139  , rows( A.rows(), A.cols(), A_trans ) );
140 
141  // Assert size matchin of B = (op(A') * M * op(A))
142  DenseLinAlgPack::Vp_V_assert_sizes(
143  B->cols(), cols( A.rows(), A.cols(), A_trans ) );
144 
145  //
146  // B = a * op(A') * M * op(A)
147  //
148  // = a * op(A') * M^(1/2) * M^(1/2) * op(A)
149  //
150  // = a * E * E'
151  //
152  // E = M^(1/2) * op(A)
153  //
154  // [ . ] [ . ]
155  // [ sqrt(M(j(1))) ] [ op(A)(j(1),:) ]
156  // [ . ] [ . ]
157  // = [ sqrt(M(j(i)) ] [ op(A)(j(i),:) ]
158  // [ . ] [ . ]
159  // [ sqrt(M(j(nz)) ] [ op(A)(j(nz),:) ]
160  // [ . ] [ . ]
161  //
162  //
163  // [ . ]
164  // [ d(j(1))' ]
165  // [ . ]
166  // = [ d(j(i))' ]
167  // [ . ]
168  // [ d(j(1))' ]
169  // [ . ]
170  //
171  // where: d(j(i)) = sqrt(M(j(i)) * op(A')(:,j(i)) <: R^m
172  // = sqrt(M(j(i)) * op(A') * e(j(i)) <: R^m
173  //
174  // Above M^(1/2) only has nz nonzero elements sqrt(M(j(i)), i = 1..nz and only
175  // the corresponding rows of op(A)(j(i),:), i = 1..nz are shown. A may in fact
176  // dense matrix but the columns are extracted through op(A)*eta(j(i)), i=1..nz.
177  //
178  // The above product B = a * E * E' is a set of nz rank-1 updates and can be written
179  // in the form:
180  //
181  // B = sum( a * d(j(i)) * d(j(i))', i = 1..nz )
182  //
183  // Since it is more efficient to perform several rank-1 updates at a time we will
184  // perform them in blocks.
185  //
186  // B = B + D(k) * D(k)', k = 1 .. num_blocks
187  //
188  // where:
189  // num_blocks = nz / num_updates_at_once + 1 (integer division)
190  // D(k) = [ d(j(i1)) ... d(j(i2)) ]
191  // i1 = (k-1) * num_updates_at_once + 1
192  // i2 = i1 + num_updates_at_once - 1
193 
194  const SpVectorSlice
195  &diag = this->diag();
196 
197  const size_type
198  n = this->rows(),
199  m = cols(A.rows(),A.cols(),A_trans);
200 
201  // Get the actual number of updates to use per rank-(num_updates) update
202  const size_type
203  num_updates
204  = my_min( num_updates_at_once()
205  ? num_updates_at_once()
206  : 20 // There may be a better default value for this?
207  , diag.nz()
208  );
209 
210  // Get the number of blocks of rank-(num_updates) updates
211  size_type
212  num_blocks = diag.nz() / num_updates;
213  if( diag.nz() % num_updates > 0 )
214  num_blocks++;
215 
216  // Initialize B = b*B
217  if( b != 1.0 )
218  assign( &nonconst_tri_ele( B->gms(), B->uplo() ), 0.0 );
219 
220  // Perform the rank-(num_updates) updates
221  DMatrix D(m,num_updates);
222  for( size_type k = 1; k <= num_blocks; ++k ) {
223  const size_type
224  i1 = (k-1) * num_updates + 1,
225  i2 = my_min( diag.nz(), i1 + num_updates - 1 );
226  // Generate the colunns of D(k)
228  m_itr = diag.begin() + (i1-1);
229  for( size_type l = 1; l <= i2-i1+1; ++l, ++m_itr ) {
230  TEUCHOS_TEST_FOR_EXCEPT( !( m_itr < diag.end() ) );
231  TEUCHOS_TEST_FOR_EXCEPT( !( m_itr->value() >= 0.0 ) );
232  V_MtV( &D.col(l), A, trans_not(A_trans)
233  , eta( m_itr->index(), n, std::sqrt(m_itr->value()) )() );
234  }
235  const DMatrixSlice
236  D_update = D(1,m,1,i2-i1+1);
237 
238 
239 // // For debugging only
240 // std::ofstream ofile("MatrixSymDiagonalSparse_Error.out");
241 // assert_print_nan_inf( D_update, "D", true, &ofile );
242  // Perform the rank-(num_updates) update
243  syrk( BLAS_Cpp::no_trans, alpha, D_update, 1.0, B );
244  }
245 }
246 
247 // Overridden from MatrixConvertToSparse
248 
249 index_type
251  EExtractRegion extract_region
252  ,EElementUniqueness element_uniqueness
253  ) const
254 {
255  return diag().nz();
256 }
257 
259  EExtractRegion extract_region
260  ,EElementUniqueness element_uniqueness
261  ,const index_type len_Aval
262  ,value_type Aval[]
263  ,const index_type len_Aij
264  ,index_type Arow[]
265  ,index_type Acol[]
266  ,const index_type row_offset
267  ,const index_type col_offset
268  ) const
269 {
270  const SpVectorSlice
271  &diag = this->diag();
272 
274  (len_Aval != 0 ? len_Aval != diag.nz() : Aval != NULL)
275  ,std::invalid_argument
276  ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
278  (len_Aij != 0 ? len_Aij != diag.nz() : (Acol != NULL || Acol != NULL) )
279  ,std::invalid_argument
280  ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
281 
282  if( len_Aval > 0 ) {
284  itr;
285  FortranTypes::f_dbl_prec
286  *l_Aval;
287  for( itr = diag.begin(), l_Aval = Aval; itr != diag.end(); ++itr ) {
288  *l_Aval++ = itr->value();
289  }
290  }
291  if( len_Aij > 0 ) {
293  itr;
294  index_type
295  *l_Arow, *l_Acol;
296  for( itr = diag.begin(), l_Arow = Arow, l_Acol = Acol; itr != diag.end(); ++itr ) {
297  const index_type
298  ij = itr->index() + diag.offset();
299  *l_Arow++ = ij + row_offset;
300  *l_Acol++ = ij + col_offset;
301  }
302  }
303 }
304 
305 } // end namespace AbstractLinAlgPack
difference_type offset() const
Return the offset for the indexes (ith real index = begin()[i-1]->index() + offset()# ...
size_type dim() const
Return the number of elements in the full vector.
virtual const SpVectorSlice diag() const =0
Give access to the sparse diagonal.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
BLAS_Cpp::Uplo uplo() const
Base class for all matrices implemented in a shared memory address space.
virtual size_type cols() const
Return the number of columns in the matrix.
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Perform a rank-k update of a symmetric matrix of the form:
DVectorSlice col(size_type j)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
Create an eta vector (scaled by alpha = default 1).
void coor_extract_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness, const index_type len_Aval, value_type Aval[], const index_type len_Aij, index_type Arow[], index_type Acol[], const index_type row_offset, const index_type col_offset) const
virtual void syrk(BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, DMatrixSliceSym *sym_lhs) const
Perform a rank-k update of a dense symmetric matrix of the form:
MatrixSymDiagSparse()
The default value of num_updates_at_once == 0 is set to allow this class to determine the appropriate...
Transp trans_not(Transp _trans)
virtual size_type rows() const
Return the number of rows in the matrix.
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
size_type nz() const
Return the number of non-zero elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void Mp_StMtMtM(DMatrixSliceSym *sym_lhs, value_type alpha, EMatRhsPlaceHolder dummy_place_holder, const MatrixOpSerial &mwo_rhs, BLAS_Cpp::Transp mwo_rhs_trans, value_type beta) const
Computes the dense symmetric matrix B += a*op(A')*M*op(A).
index_type num_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness) const