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