MFEM  v3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
table.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of data types Table.
13 
14 #include <iostream>
15 #include <iomanip>
16 
17 #include "array.hpp"
18 #include "table.hpp"
19 #include "error.hpp"
20 
21 namespace mfem
22 {
23 
24 using namespace std;
25 
26 Table::Table(const Table &table)
27 {
28  size = table.size;
29  if (size >= 0)
30  {
31  const int nnz = table.I[size];
32  I = new int[size+1];
33  J = new int[nnz];
34  memcpy(I, table.I, sizeof(int)*(size+1));
35  memcpy(J, table.J, sizeof(int)*nnz);
36  }
37  else
38  {
39  I = J = NULL;
40  }
41 }
42 
43 Table::Table (int dim, int connections_per_row)
44 {
45  int i, j, sum = dim * connections_per_row;
46 
47  size = dim;
48  I = new int[size+1];
49  J = new int[sum];
50 
51  I[0] = 0;
52  for (i = 1; i <= size; i++)
53  {
54  I[i] = I[i-1] + connections_per_row;
55  for (j = I[i-1]; j < I[i]; j++) { J[j] = -1; }
56  }
57 }
58 
59 Table::Table (int nrows, int *partitioning)
60 {
61  size = nrows;
62 
63  I = new int[size+1];
64  J = new int[size];
65 
66  for (int i = 0; i < size; i++)
67  {
68  I[i] = i;
69  J[i] = partitioning[i];
70  }
71  I[size] = size;
72 }
73 
74 void Table::MakeI (int nrows)
75 {
76  SetDims (nrows, 0);
77 
78  for (int i = 0; i <= nrows; i++)
79  {
80  I[i] = 0;
81  }
82 }
83 
85 {
86  int i, j, k;
87 
88  for (k = i = 0; i < size; i++)
89  {
90  j = I[i], I[i] = k, k += j;
91  }
92 
93  J = new int[I[size]=k];
94 }
95 
96 void Table::AddConnections (int r, const int *c, int nc)
97 {
98  int *jp = J+I[r];
99 
100  for (int i = 0; i < nc; i++)
101  {
102  jp[i] = c[i];
103  }
104  I[r] += nc;
105 }
106 
108 {
109  for (int i = size; i > 0; i--)
110  {
111  I[i] = I[i-1];
112  }
113  I[0] = 0;
114 }
115 
116 void Table::SetSize(int dim, int connections_per_row)
117 {
118  SetDims (dim, dim * connections_per_row);
119 
120  if (size > 0)
121  {
122  I[0] = 0;
123  for (int i = 0, j = 0; i < size; i++)
124  {
125  int end = I[i] + connections_per_row;
126  I[i+1] = end;
127  for ( ; j < end; j++) { J[j] = -1; }
128  }
129  }
130 }
131 
132 void Table::SetDims(int rows, int nnz)
133 {
134  int j;
135 
136  j = (I) ? (I[size]) : (0);
137  if (size != rows)
138  {
139  size = rows;
140  if (I) { delete [] I; }
141  I = (rows >= 0) ? (new int[rows+1]) : (NULL);
142  }
143 
144  if (j != nnz)
145  {
146  if (J) { delete [] J; }
147  J = (nnz > 0) ? (new int[nnz]) : (NULL);
148  }
149 
150  if (size >= 0)
151  {
152  I[0] = 0;
153  I[size] = nnz;
154  }
155 }
156 
157 int Table::operator() (int i, int j) const
158 {
159  if ( i>=size || i<0 )
160  {
161  return -1;
162  }
163 
164  int k, end = I[i+1];
165  for (k = I[i]; k < end; k++)
166  {
167  if (J[k] == j)
168  {
169  return k;
170  }
171  else if (J[k] == -1)
172  {
173  return -1;
174  }
175  }
176  return -1;
177 }
178 
179 void Table::GetRow(int i, Array<int> &row) const
180 {
181  MFEM_ASSERT(i >= 0 && i < size, "Row index " << i << " is out of range [0,"
182  << size << ')');
183 
184  row.SetSize(RowSize(i));
185  row.Assign(GetRow(i));
186 }
187 
189 {
190  for (int r = 0; r < size; r++)
191  {
192  std::sort(J + I[r], J + I[r+1]);
193  }
194 }
195 
196 void Table::SetIJ(int *newI, int *newJ, int newsize)
197 {
198  delete [] I;
199  delete [] J;
200  I = newI;
201  J = newJ;
202  if (newsize >= 0)
203  {
204  size = newsize;
205  }
206 }
207 
208 int Table::Push(int i, int j)
209 {
210  MFEM_ASSERT( i >=0 && i<size, "Index out of bounds. i = "<<i);
211 
212  for (int k = I[i], end = I[i+1]; k < end; k++)
213  if (J[k] == j)
214  {
215  return k;
216  }
217  else if (J[k] == -1)
218  {
219  J[k] = j;
220  return k;
221  }
222 
223  MFEM_ABORT("Reached end of loop unexpectedly: (i,j) = (" << i << ", " << j
224  << ")");
225 
226  return -1;
227 }
228 
230 {
231  int i, j, end, sum = 0, n = 0, newI = 0;
232 
233  for (i=0; i<I[size]; i++)
234  if (J[i] != -1)
235  {
236  sum++;
237  }
238 
239  if (sum != I[size])
240  {
241  int *NewJ = new int[sum];
242 
243  for (i=0; i<size; i++)
244  {
245  end = I[i+1];
246  for (j=I[i]; j<end; j++)
247  {
248  if (J[j] == -1) { break; }
249  NewJ[ n++ ] = J[j];
250  }
251  I[i] = newI;
252  newI = n;
253  }
254  I[size] = sum;
255 
256  delete [] J;
257 
258  J = NewJ;
259 
260  MFEM_ASSERT(sum == n, "sum = " << sum << ", n = " << n);
261  }
262 }
263 
264 void Table::MakeFromList(int nrows, const Array<Connection> &list)
265 {
266  Clear();
267 
268  size = nrows;
269  int nnz = list.Size();
270 
271  I = new int[size+1];
272  J = new int[nnz];
273 
274  for (int i = 0, k = 0; i <= size; i++)
275  {
276  I[i] = k;
277  while (k < nnz && list[k].from == i)
278  {
279  J[k] = list[k].to;
280  k++;
281  }
282  }
283 }
284 
285 int Table::Width() const
286 {
287  int width = -1, nnz = (size >= 0) ? I[size] : 0;
288  for (int k = 0; k < nnz; k++)
289  {
290  if (J[k] > width) { width = J[k]; }
291  }
292  return width + 1;
293 }
294 
295 void Table::Print(std::ostream & out, int width) const
296 {
297  int i, j;
298 
299  for (i = 0; i < size; i++)
300  {
301  out << "[row " << i << "]\n";
302  for (j = I[i]; j < I[i+1]; j++)
303  {
304  out << setw(5) << J[j];
305  if ( !((j+1-I[i]) % width) )
306  {
307  out << '\n';
308  }
309  }
310  if ((j-I[i]) % width)
311  {
312  out << '\n';
313  }
314  }
315 }
316 
317 void Table::PrintMatlab(std::ostream & out) const
318 {
319  int i, j;
320 
321  for (i = 0; i < size; i++)
322  for (j = I[i]; j < I[i+1]; j++)
323  {
324  out << i << " " << J[j] << " 1. \n";
325  }
326 
327  out << flush;
328 }
329 
330 void Table::Save(std::ostream & out) const
331 {
332  int i;
333 
334  out << size << '\n';
335 
336  for (i = 0; i <= size; i++)
337  {
338  out << I[i] << '\n';
339  }
340  for (i = 0; i < I[size]; i++)
341  {
342  out << J[i] << '\n';
343  }
344 }
345 
347 {
348  delete [] I;
349  delete [] J;
350  size = -1;
351  I = J = NULL;
352 }
353 
354 void Table::Copy(Table & copy) const
355 {
356  if (size >= 0)
357  {
358  int * i_copy = new int[size+1];
359  int * j_copy = new int[I[size]];
360 
361  memcpy(i_copy, I, sizeof(int)*(size+1));
362  memcpy(j_copy, J, sizeof(int)*I[size]);
363 
364  copy.SetIJ(i_copy, j_copy, size);
365  }
366  else
367  {
368  copy.Clear();
369  }
370 }
371 
372 void Table::Swap(Table & other)
373 {
374  mfem::Swap(size, other.size);
375  mfem::Swap(I, other.I);
376  mfem::Swap(J, other.J);
377 }
378 
379 long Table::MemoryUsage() const
380 {
381  if (size < 0 || I == NULL) { return 0; }
382  return (size+1 + I[size]) * sizeof(int);
383 }
384 
386 {
387  if (I) { delete [] I; }
388  if (J) { delete [] J; }
389 }
390 
391 void Transpose (const Table &A, Table &At, int _ncols_A)
392 {
393  const int *i_A = A.GetI();
394  const int *j_A = A.GetJ();
395  const int nrows_A = A.Size();
396  const int ncols_A = (_ncols_A < 0) ? A.Width() : _ncols_A;
397  const int nnz_A = i_A[nrows_A];
398 
399  At.SetDims (ncols_A, nnz_A);
400 
401  int *i_At = At.GetI();
402  int *j_At = At.GetJ();
403 
404  for (int i = 0; i <= ncols_A; i++)
405  {
406  i_At[i] = 0;
407  }
408  for (int i = 0; i < nnz_A; i++)
409  {
410  i_At[j_A[i]+1]++;
411  }
412  for (int i = 1; i < ncols_A; i++)
413  {
414  i_At[i+1] += i_At[i];
415  }
416 
417  for (int i = 0; i < nrows_A; i++)
418  for (int j = i_A[i]; j < i_A[i+1]; j++)
419  {
420  j_At[i_At[j_A[j]]++] = i;
421  }
422  for (int i = ncols_A; i > 0; i--)
423  {
424  i_At[i] = i_At[i-1];
425  }
426  i_At[0] = 0;
427 }
428 
429 Table * Transpose(const Table &A)
430 {
431  Table * At = new Table;
432  Transpose(A, *At);
433  return At;
434 }
435 
436 void Transpose(const Array<int> &A, Table &At, int _ncols_A)
437 {
438  At.MakeI((_ncols_A < 0) ? (A.Max() + 1) : _ncols_A);
439  for (int i = 0; i < A.Size(); i++)
440  {
441  At.AddAColumnInRow(A[i]);
442  }
443  At.MakeJ();
444  for (int i = 0; i < A.Size(); i++)
445  {
446  At.AddConnection(A[i], i);
447  }
448  At.ShiftUpI();
449 }
450 
451 void Mult (const Table &A, const Table &B, Table &C)
452 {
453  int i, j, k, l, m;
454  const int *i_A = A.GetI();
455  const int *j_A = A.GetJ();
456  const int *i_B = B.GetI();
457  const int *j_B = B.GetJ();
458  const int nrows_A = A.Size();
459  const int nrows_B = B.Size();
460  const int ncols_A = A.Width();
461  const int ncols_B = B.Width();
462 
463  MFEM_VERIFY( ncols_A <= nrows_B, "Table size mismatch: ncols_A = " << ncols_A
464  << ", nrows_B = " << nrows_B);
465 
466  Array<int> B_marker (ncols_B);
467 
468  for (i = 0; i < ncols_B; i++)
469  {
470  B_marker[i] = -1;
471  }
472 
473  int counter = 0;
474  for (i = 0; i < nrows_A; i++)
475  {
476  for (j = i_A[i]; j < i_A[i+1]; j++)
477  {
478  k = j_A[j];
479  for (l = i_B[k]; l < i_B[k+1]; l++)
480  {
481  m = j_B[l];
482  if (B_marker[m] != i)
483  {
484  B_marker[m] = i;
485  counter++;
486  }
487  }
488  }
489  }
490 
491  C.SetDims (nrows_A, counter);
492 
493  for (i = 0; i < ncols_B; i++)
494  {
495  B_marker[i] = -1;
496  }
497 
498  int *i_C = C.GetI();
499  int *j_C = C.GetJ();
500  counter = 0;
501  for (i = 0; i < nrows_A; i++)
502  {
503  i_C[i] = counter;
504  for (j = i_A[i]; j < i_A[i+1]; j++)
505  {
506  k = j_A[j];
507  for (l = i_B[k]; l < i_B[k+1]; l++)
508  {
509  m = j_B[l];
510  if (B_marker[m] != i)
511  {
512  B_marker[m] = i;
513  j_C[counter] = m;
514  counter++;
515  }
516  }
517  }
518  }
519 }
520 
521 
522 Table * Mult (const Table &A, const Table &B)
523 {
524  Table * C = new Table;
525  Mult(A,B,*C);
526  return C;
527 }
528 
529 STable::STable (int dim, int connections_per_row) :
530  Table(dim, connections_per_row)
531 {}
532 
533 int STable::operator() (int i, int j) const
534 {
535  if (i < j)
536  {
537  return Table::operator()(i,j);
538  }
539  else
540  {
541  return Table::operator()(j,i);
542  }
543 }
544 
545 int STable::Push( int i, int j )
546 {
547  if (i < j)
548  {
549  return Table::Push(i, j);
550  }
551  else
552  {
553  return Table::Push(j, i);
554  }
555 }
556 
557 
559 {
560  Rows = new Node*[nrows];
561  for (int i = 0; i < nrows; i++)
562  {
563  Rows[i] = NULL;
564  }
565  NumRows = nrows;
566  NumEntries = 0;
567 }
568 
569 int DSTable::Push_(int r, int c)
570 {
571  MFEM_ASSERT(r >= 0 && r < NumRows,
572  "Row out of bounds: r = " << r << ", NumRows = " << NumRows);
573  Node *n;
574  for (n = Rows[r]; n != NULL; n = n->Prev)
575  {
576  if (n->Column == c)
577  {
578  return (n->Index);
579  }
580  }
581 #ifdef MFEM_USE_MEMALLOC
582  n = NodesMem.Alloc ();
583 #else
584  n = new Node;
585 #endif
586  n->Column = c;
587  n->Index = NumEntries;
588  n->Prev = Rows[r];
589  Rows[r] = n;
590  return (NumEntries++);
591 }
592 
593 int DSTable::Index(int r, int c) const
594 {
595  MFEM_ASSERT( r>=0, "Row index must be non-negative, not "<<r);
596  if (r >= NumRows)
597  {
598  return (-1);
599  }
600  for (Node *n = Rows[r]; n != NULL; n = n->Prev)
601  {
602  if (n->Column == c)
603  {
604  return (n->Index);
605  }
606  }
607  return (-1);
608 }
609 
611 {
612 #ifdef MFEM_USE_MEMALLOC
613  // NodesMem.Clear(); // this is done implicitly
614 #else
615  for (int i = 0; i < NumRows; i++)
616  {
617  Node *na, *nb = Rows[i];
618  while (nb != NULL)
619  {
620  na = nb;
621  nb = nb->Prev;
622  delete na;
623  }
624  }
625 #endif
626  delete [] Rows;
627 }
628 
629 }
Elem * Alloc()
Definition: mem_alloc.hpp:124
int Size() const
Logical size of the array.
Definition: array.hpp:109
int Push(int i, int j)
Definition: table.cpp:208
int * GetJ()
Definition: table.hpp:108
void SetSize(int dim, int connections_per_row)
Set the size and the number of connections for the table.
Definition: table.cpp:116
int operator()(int i, int j) const
Definition: table.cpp:157
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:74
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition: table.cpp:188
void Swap(Table &other)
Definition: table.cpp:372
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:451
void SetDims(int rows, int nnz)
Definition: table.cpp:132
int * J
Definition: table.hpp:49
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
int Push(int i, int j)
Definition: table.cpp:545
void Save(std::ostream &out) const
Definition: table.cpp:330
void AddConnections(int r, const int *c, int nc)
Definition: table.cpp:96
int Width() const
Returns the number of TYPE II elements (after Finalize() is called).
Definition: table.cpp:285
void Print(std::ostream &out=std::cout, int width=4) const
Prints the table to stream out.
Definition: table.cpp:295
int operator()(int i, int j) const
Definition: table.cpp:533
int dim
Definition: ex3.cpp:47
void MakeFromList(int nrows, const Array< Connection > &list)
Definition: table.cpp:264
void Clear()
Definition: table.cpp:346
void AddConnection(int r, int c)
Definition: table.hpp:74
int * I
Definition: table.hpp:49
T Max() const
Definition: array.cpp:90
long MemoryUsage() const
Definition: table.cpp:379
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:391
void Assign(const T *)
Copy data from a pointer. Size() elements are copied.
Definition: array.hpp:518
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
void PrintMatlab(std::ostream &out) const
Definition: table.cpp:317
STable(int dim, int connections_per_row=3)
Creates table with fixed number of connections.
Definition: table.cpp:529
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:322
void Finalize()
Definition: table.cpp:229
int size
size is the number of TYPE I elements.
Definition: table.hpp:44
void AddAColumnInRow(int r)
Definition: table.hpp:71
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
void ShiftUpI()
Definition: table.cpp:107
~Table()
Destroys Table.
Definition: table.cpp:385
void SetIJ(int *newI, int *newJ, int newsize=-1)
Definition: table.cpp:196
void MakeJ()
Definition: table.cpp:84
int * GetI()
Definition: table.hpp:107
void Copy(Table &copy) const
Definition: table.cpp:354
DSTable(int nrows)
Definition: table.cpp:558
Table()
Creates an empty table.
Definition: table.hpp:53