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_MatrixConvertToSparseEncap.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 "AbstractLinAlgPack_MatrixConvertToSparseEncap.hpp"
45 #include "AbstractLinAlgPack_MatrixExtractSparseElements.hpp"
46 #include "AbstractLinAlgPack_VectorSpace.hpp"
47 #include "DenseLinAlgPack_IVector.hpp"
48 #include "Teuchos_Assert.hpp"
49 
50 namespace AbstractLinAlgPack {
51 
52 // Constructors/initializers
53 
55  :row_rng_(Range1D::Invalid)
56  ,col_rng_(Range1D::Invalid)
57  ,mese_trans_(BLAS_Cpp::no_trans)
58  ,alpha_(0.0)
59  ,nz_full_(0)
60 {}
61 
63  const mese_ptr_t &mese
64  ,const i_vector_ptr_t &inv_row_perm
65  ,const Range1D &row_rng
66  ,const i_vector_ptr_t &inv_col_perm
67  ,const Range1D &col_rng
68  ,const BLAS_Cpp::Transp mese_trans
69  ,const value_type alpha
70  )
71 {
72  this->initialize(mese,inv_row_perm,row_rng,inv_col_perm,col_rng,mese_trans,alpha);
73 }
74 
76  const mese_ptr_t &mese
77  ,const i_vector_ptr_t &inv_row_perm
78  ,const Range1D &row_rng_in
79  ,const i_vector_ptr_t &inv_col_perm
80  ,const Range1D &col_rng_in
81  ,const BLAS_Cpp::Transp mese_trans
82  ,const value_type alpha
83  )
84 {
85  const size_type mese_rows = mese->rows(), mese_cols = mese->cols();
86  const Range1D row_rng = RangePack::full_range(row_rng_in,1,mese_rows);
87  const Range1D col_rng = RangePack::full_range(col_rng_in,1,mese_cols);
88 #ifdef TEUCHOS_DEBUG
89  const char msg_head[] = "MatrixConvertToSparseEncap::initialize(...): Error!";
90  TEUCHOS_TEST_FOR_EXCEPTION( mese.get() == NULL, std::logic_error, msg_head );
91  TEUCHOS_TEST_FOR_EXCEPTION( inv_row_perm.get() != NULL && inv_row_perm->size() != mese_rows, std::logic_error, msg_head );
92  TEUCHOS_TEST_FOR_EXCEPTION( row_rng.ubound() > mese_rows, std::logic_error, msg_head );
93  TEUCHOS_TEST_FOR_EXCEPTION( inv_col_perm.get() != NULL && inv_col_perm->size() != mese_cols, std::logic_error, msg_head );
94  TEUCHOS_TEST_FOR_EXCEPTION( col_rng.ubound() > mese->cols(), std::logic_error, msg_head );
95 #endif
96  mese_ = mese;
97  mese_trans_ = mese_trans;
98  alpha_ = alpha;
99  inv_row_perm_ = inv_row_perm;
100  inv_col_perm_ = inv_col_perm;
101  row_rng_ = row_rng;
102  col_rng_ = col_rng;
104  space_cols_ = ( mese_trans_ == BLAS_Cpp::no_trans
105  ? mese_->space_cols().sub_space(row_rng_)
106  : mese_->space_rows().sub_space(col_rng_) );
107  space_rows_ = ( mese_trans_ == BLAS_Cpp::no_trans
108  ? mese_->space_rows().sub_space(col_rng_)
109  : mese_->space_cols().sub_space(row_rng_) );
110 }
111 
113 {
114  namespace mmp = MemMngPack;
115  mese_ = Teuchos::null;
116  inv_row_perm_ = Teuchos::null;
117  row_rng_ = Range1D::Invalid;
118  inv_col_perm_ = Teuchos::null;
119  col_rng_ = Range1D::Invalid;
120  mese_trans_ = BLAS_Cpp::no_trans;
121  alpha_ = 0.0;
122  nz_full_ = 0;
123 }
124 
125 // Overridden from MatrixBase
126 
128 {
129  return *space_cols_;
130 }
131 
133 {
134  return *space_rows_;
135 }
136 
138 {
139  return mese_trans_ == BLAS_Cpp::no_trans ? row_rng_.size() : col_rng_.size();
140 }
141 
143 {
144  return mese_trans_ == BLAS_Cpp::no_trans ? col_rng_.size() : row_rng_.size();
145 }
146 
148 {
149  return nz_full_;
150 }
151 
152 // Overridden from MatrixConvertToSparse
153 
155  EExtractRegion extract_region_in
156  ,EElementUniqueness element_uniqueness
157  ) const
158 {
159  index_type dl = 0, du = 0;
160  EExtractRegion extract_region = extract_region_in;
161  if( mese_trans_ == BLAS_Cpp::trans )
162  extract_region
163  = ( extract_region_in == EXTRACT_FULL_MATRIX
165  : ( extract_region_in == EXTRACT_UPPER_TRIANGULAR
168  switch(extract_region) {
169  case EXTRACT_FULL_MATRIX:
170  dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
171  du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
172  break;
174  dl = -(index_type)row_rng_.lbound() + (index_type)col_rng_.lbound();
175  du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
176  break;
178  dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
179  du = +(index_type)col_rng_.lbound() - (index_type)row_rng_.lbound();
180  break;
181  default:
183  }
184  const index_type
185  *inv_row_perm = inv_row_perm_.get() ? &(*inv_row_perm_)(1) : NULL,
186  *inv_col_perm = inv_col_perm_.get() ? &(*inv_col_perm_)(1) : NULL;
187  return mese_->count_nonzeros(
188  element_uniqueness
189  ,inv_row_perm
190  ,inv_col_perm
191  ,row_rng_
192  ,col_rng_
193  ,dl
194  ,du
195  );
196 }
197 
199  EExtractRegion extract_region
200  ,EElementUniqueness element_uniqueness
201  ,const index_type len_Aval
202  ,value_type Aval[]
203  ,const index_type len_Aij
204  ,index_type Arow[]
205  ,index_type Acol[]
206  ,const index_type row_offset
207  ,const index_type col_offset
208  ) const
209 {
210  index_type dl = 0, du = 0; // This may not be correct!
211  switch(extract_region) {
212  case EXTRACT_FULL_MATRIX:
213  dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
214  du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
215  break;
217  dl = -(index_type)row_rng_.lbound() + (index_type)col_rng_.lbound();
218  du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
219  break;
221  dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
222  du = +(index_type)col_rng_.lbound() - (index_type)row_rng_.lbound();
223  break;
224  default:
226  }
227  const index_type
228  *inv_row_perm = inv_row_perm_.get() ? &(*inv_row_perm_)(1) : NULL,
229  *inv_col_perm = inv_col_perm_.get() ? &(*inv_col_perm_)(1) : NULL;
230  mese_->coor_extract_nonzeros(
231  element_uniqueness
232  ,inv_row_perm
233  ,inv_col_perm
234  ,row_rng_
235  ,col_rng_
236  ,dl
237  ,du
238  ,alpha_
239  ,len_Aval
240  ,Aval
241  ,len_Aij
242  ,mese_trans_ == BLAS_Cpp::no_trans ? Arow : Acol
243  ,mese_trans_ == BLAS_Cpp::no_trans ? Acol : Arow
244  ,mese_trans_ == BLAS_Cpp::no_trans ? row_offset : col_offset
245  ,mese_trans_ == BLAS_Cpp::no_trans ? col_offset : row_offset
246  );
247 }
248 
249 } // end namespace AbstractLinAlgPack
const i_vector_ptr_t & inv_row_perm() const
index_type num_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness) const
const i_vector_ptr_t & inv_col_perm() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index ubound() const
T * get() const
Abstract interface for objects that represent a space for mutable coordinate vectors.
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
Entries with duplicate row and column indexes are allowed with the understanding that the values are ...
Extract the nonzero elements from the lower triangular region including the diagonal.
static const Range1D Invalid
void initialize(const mese_ptr_t &mese, const i_vector_ptr_t &inv_row_perm, const Range1D &row_rng, const i_vector_ptr_t &inv_col_perm, const Range1D &col_rng, const BLAS_Cpp::Transp mese_trans, const value_type alpha=1.0)
Initialize a permuted view of a sparse matrix.
Transp
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Extract the nonzero elements from the upper triangular region including the diagonal.