Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
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...
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
void get_node_ids(int elem_id, int *node_ids)
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumMyElements() const
Number of elements on the calling processor.
matrix_data(int num_rows, int *rowlengths, int blocksize=1)
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.
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
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)
long long * MyGlobalElements64() const
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