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