Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_matrix_data.cpp
Go to the documentation of this file.
1
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42
43#include <Epetra_matrix_data.h>
44#include <Epetra_Map.h>
45#include <Epetra_CrsMatrix.h>
46#include <Epetra_Util.h>
47
48namespace epetra_test {
49
51 int* rowlengths_in,
52 int blocksize_in)
53 : numrows_(num_rows),
54 numcols_(0),
55 rows_(0),
56 rowlengths_(0),
57 blocksize_(blocksize_in),
58 colindices_(0),
59 coefs_(0)
60{
61 if (numrows_ > 0) {
62 rows_ = new int[numrows_];
63 rowlengths_ = new int[numrows_];
64 colindices_ = new int*[numrows_];
65 coefs_ = new double*[numrows_];
66 int dim = blocksize_in*blocksize_in;
67 for(int i=0; i<numrows_; ++i) {
68 rows_[i] = i;
69 rowlengths_[i] = rowlengths_in[i];
70 colindices_[i] = new int[rowlengths_[i]];
71 coefs_[i] = new double[rowlengths_[i]*dim];
72
73 for(int j=0; j<rowlengths_[i]; ++j) {
74 colindices_[i][j] = 0;
75 for(int k=0; k<dim; ++k) coefs_[i][j*dim+k] = 0.0;
76 }
77 }
78 }
79}
80
82 int num_cols,
83 int num_off_diagonals,
84 int blocksize_in)
85 : numrows_(num_rows),
86 numcols_(num_cols),
87 rows_(0),
88 rowlengths_(0),
89 blocksize_(blocksize_in),
90 colindices_(0),
91 coefs_(0)
92{
93 if (numrows_ > 0) {
94 rows_ = new int[numrows_];
95 rowlengths_ = new int[numrows_];
96 colindices_ = new int*[numrows_];
97 coefs_ = new double*[numrows_];
98
99 int max_row_length = 1+num_off_diagonals*2;
100
101 for(int i=0; i<numrows_; ++i) {
102 rows_[i] = i;
103 if (i >= num_off_diagonals && numrows_-i > num_off_diagonals) {
104 rowlengths_[i] = max_row_length;
105 }
106 else {
107 if (i<num_off_diagonals) {
108 rowlengths_[i] = 1+max_row_length/2+i;
109 }
110 else {
111 rowlengths_[i] = 1+max_row_length/2+numrows_-i-1;
112 }
113 }
114 colindices_[i] = new int[rowlengths_[i]];
115 int dim = blocksize_in*blocksize_in;
116 coefs_[i] = new double[rowlengths_[i]*dim];
117
118 int first_col = i - max_row_length/2;
119 if (first_col < 0) first_col = 0;
120
121 for(int j=0; j<rowlengths_[i]; ++j) {
122 colindices_[i][j] = first_col+j;
123 for(int k=0; k<dim; ++k) {
124 coefs_[i][j*dim+k] = 1.0;
125 }
126 }
127 }
128 }
129}
130
131static const int nodes_per_elem = 4;
132
133void get_node_ids(int elem_id, int* node_ids)
134{
135 int first_node = 2*elem_id;
136 for(int i=0; i<nodes_per_elem; ++i) node_ids[i] = first_node+i;
137}
138
139matrix_data::matrix_data(int num_quad_elements,
140 int num_dof_per_node,
141 bool make_numerically_nonsymmetric)
142 : numrows_(0),
143 numcols_(0),
144 rows_(0),
145 rowlengths_(0),
146 blocksize_(num_dof_per_node),
147 colindices_(0),
148 coefs_(0)
149{
150 //Set up matrix-data representing a simple finite-element
151 //mesh containing 2-D quad elements
152 //
153 // *-----*-----*-----*
154 // 0| 2| 4| 6|
155 // | 0 | 1 | ne-1|
156 // | | | |
157 // *-----*-----*-----*
158 // 1 3 5 7
159 //
160 //In the above drawing, 'ne' means num-elems. node-numbers are to the
161 //lower-left of each node (*).
162
163 numrows_ = num_quad_elements*2+2;
164
165 if (numrows_ > 0) {
166 rows_ = new int[numrows_];
167 rowlengths_ = new int[numrows_];
168 colindices_ = new int*[numrows_];
169 coefs_ = new double*[numrows_];
170
171 int i, j, k;
172 for(i=0; i<numrows_; ++i) {
173 rows_[i] = i;
174 rowlengths_[i] = 0;
175 }
176
177 int* nodes = new int[nodes_per_elem];
178 for(i=0; i<num_quad_elements; ++i) {
179 get_node_ids(i, nodes);
180
181 for(j=0; j<nodes_per_elem; ++j) {
182 int node_j = nodes[j];
183 for(k=0; k<nodes_per_elem; ++k) {
184 int insertPoint = -1;
185 int alloclen = rowlengths_[node_j];
186 int offset = Epetra_Util_binary_search(nodes[k], colindices_[node_j],
187 rowlengths_[node_j], insertPoint);
188 if (offset < 0) {
189 Epetra_Util_insert(nodes[k], insertPoint,
190 colindices_[node_j], rowlengths_[node_j],
191 alloclen);
192 }
193 }
194 }
195 }
196
197 int dim = blocksize_*blocksize_;
198 for(i=0; i<numrows_; ++i) {
199 int len = rowlengths_[i]*dim;
200 coefs_[i] = new double[len];
201 for(j=0; j<len; ++j) {
202 if (make_numerically_nonsymmetric) {
203 coefs_[i][j] = 1.0*(j+1);
204 }
205 else {
206 coefs_[i][j] = 1.0;
207 }
208 }
209 }
210 }
211}
212
214{
215 for(int i=0; i<numrows_; ++i) {
216 delete [] colindices_[i];
217 delete [] coefs_[i];
218 }
219
220 delete [] colindices_; colindices_ = 0;
221 delete [] coefs_; coefs_ = 0;
222 delete [] rowlengths_; rowlengths_ = 0;
223 delete [] rows_; rows_ = 0;
224 numrows_ = 0;
225}
226
227double* matrix_data::coefs(int row, int col)
228{
229 int insertPoint = -1;
230 int row_idx = Epetra_Util_binary_search(row, rows_, numrows_,
231 insertPoint);
232 if (row_idx < 0) {
233 std::cerr << "ERROR, row " << row
234 << " not found in matrix_data"<< std::endl;
235 return 0;
236 }
237
238 int col_idx = Epetra_Util_binary_search(col, colindices_[row_idx],
239 rowlengths_[row_idx], insertPoint);
240 if (col_idx < 0) {
241 std::cerr << "ERROR, col " << col
242 << " not found in matrix_data"<< std::endl;
243 return 0;
244 }
245
246 int dim = blocksize_*blocksize_;
247 return( &(coefs_[row_idx][col_idx*dim]) );
248}
249
250/* Not used
251void matrix_data::copy_local_data_to_matrix(Epetra_CrsMatrix& A)
252{
253 const Epetra_Map& rowmap = A.RowMap();
254
255 for(int i=0; i<numrows_; ++i) {
256 int row = rows_[i];
257 if (rowmap.MyGID(row)) {
258 int err = A.ReplaceGlobalValues(row, rowlengths_[i],
259 coefs_[i], colindices_[i]);
260 if (err < 0) {
261 err = A.InsertGlobalValues(row, rowlengths_[i],
262 coefs_[i], colindices_[i]);
263 }
264 }
265 }
266}
267*/
268
270{
271 const Epetra_Map& map = A.RowMap();
272 int numMyRows = map.NumMyElements();
273 Epetra_Util util;
274
275 if(map.GlobalIndicesInt()) {
276#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
277 int* myRows_int = map.MyGlobalElements();
278 for(int i=0; i<numMyRows; ++i) {
279 int row = myRows_int[i];
280 int rowLen = A.NumGlobalEntries(row);
281 if (rowLen != rowlengths_[row]) {
282 return(false);
283 }
284
285 int* indices = new int[rowLen];
286 double* values = new double[rowLen];
287 A.ExtractGlobalRowCopy(row, rowLen, rowLen, values, indices);
288
289 util.Sort(true, rowLen, indices, 1, &values, 0, 0, 0, 0);
290
291 bool same = true;
292 int* this_indices = colindices_[row];
293 double* this_coefs = coefs_[row];
294
295 for(int j=0; j<rowLen; ++j) {
296 if (indices[j] != this_indices[j]) {
297 same = false; break;
298 }
299 if (values[j] != this_coefs[j]) {
300 same = false; break;
301 }
302 }
303
304 delete [] indices;
305 delete [] values;
306
307 if (!same) return(false);
308 }
309
310#else
311 throw "matrix_data::compare_local_data: global index int but no API for it.";
312#endif
313 }
314 else if(map.GlobalIndicesLongLong()) {
315#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
316 long long* myRows_LL = map.MyGlobalElements64();
317 for(int i=0; i<numMyRows; ++i) {
318 long long row = myRows_LL[i];
319 int rowLen = A.NumGlobalEntries(row);
320 if (rowLen != rowlengths_[row]) {
321 return(false);
322 }
323
324 long long* indices = new long long[rowLen];
325 double* values = new double[rowLen];
326 A.ExtractGlobalRowCopy(row, rowLen, rowLen, values, indices);
327
328 util.Sort(true, rowLen, indices, 1, &values, 0, 0, 0, 0);
329
330 bool same = true;
331 int* this_indices = colindices_[row];
332 double* this_coefs = coefs_[row];
333
334 for(int j=0; j<rowLen; ++j) {
335 if (indices[j] != this_indices[j]) {
336 same = false; break;
337 }
338 if (values[j] != this_coefs[j]) {
339 same = false; break;
340 }
341 }
342
343 delete [] indices;
344 delete [] values;
345
346 if (!same) return(false);
347 }
348
349#else
350 throw "matrix_data::compare_local_data: global index long long but no API for it.";
351#endif
352 }
353 else {
354 assert(false);
355 throw "matrix_data::compare_local_data: global index type of map unknown.";
356 }
357
358 return(true);
359}
360
361}//namespace epetra_test
362
int Epetra_Util_binary_search(T item, const T *list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
int Epetra_Util_insert(T item, int offset, T *&list, int &usedLength, int &allocatedLength, int allocChunkSize=32)
Function to insert an item in a list, at a specified offset.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
int NumMyElements() const
Number of elements on the calling processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_Util: The Epetra Util Wrapper Class.
Definition Epetra_Util.h:79
static void EPETRA_LIB_DLL_EXPORT Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_Util Sort Routine (Shell sort)
bool compare_local_data(const Epetra_CrsMatrix &A)
The portion of this matrix_data object's data that corresponds to the locally-owned rows of A,...
matrix_data(int num_rows, int *rowlengths, int blocksize=1)
static const int nodes_per_elem
void get_node_ids(int elem_id, int *node_ids)