Epetra Package Browser (Single Doxygen Collection)  Development
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 
48 namespace 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 
131 static const int nodes_per_elem = 4;
132 
133 void 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 
139 matrix_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 
227 double* 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
251 void 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 
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
bool compare_local_data(const Epetra_CrsMatrix &A)
The portion of this matrix_data object&#39;s data that corresponds to the locally-owned rows of A...
long long * MyGlobalElements64() const
void get_node_ids(int elem_id, int *node_ids)
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
matrix_data(int num_rows, int *rowlengths, int blocksize=1)
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
int NumMyElements() const
Number of elements on the calling processor.
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
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.
static const int nodes_per_elem
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_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
static void 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 GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
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.
Definition: Epetra_Util.h:398