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