gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
linearSystemCSR.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <stdlib.h>
7 #include <stdio.h>
8 #include <string.h>
9 #include <complex>
10 #include "GmshConfig.h"
11 #include "GmshMessage.h"
12 #include "linearSystemCSR.h"
13 #include "OS.h"
14 
15 #define SWAP(a, b) \
16  temp = (a); \
17  (a) = (b); \
18  (b) = temp;
19 #define SWAPI(a, b) \
20  tempi = (a); \
21  (a) = (b); \
22  (b) = tempi;
23 
24 static void *CSRMalloc(size_t size)
25 {
26  void *ptr;
27  if(!size) return (nullptr);
28  ptr = malloc(size);
29  return (ptr);
30 }
31 
32 static void *CSRRealloc(void *ptr, size_t size) { return realloc(ptr, size); }
33 
34 static void CSRList_Realloc(CSRList_T *liste, int n)
35 {
36  char *temp;
37  if(n <= 0) return;
38  if(liste->array == nullptr) {
39  liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
40  liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
41  }
42  else {
43  if(n > liste->nmax) {
44  liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
45  temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
46  liste->array = temp;
47  }
48  }
49 }
50 
51 static void CSRList_Resize_strict(CSRList_T *liste, int n)
52 {
53  liste->array = (char *)CSRRealloc(liste->array, n * liste->size);
54  liste->n = n;
55  liste->nmax = n;
56 }
57 
58 static CSRList_T *CSRList_Create(int n, int incr, int size)
59 {
60  CSRList_T *liste;
61 
62  if(n <= 0) n = 1;
63  if(incr <= 0) incr = 1;
64 
65  liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
66 
67  liste->nmax = 0;
68  liste->incr = incr;
69  liste->size = size;
70  liste->n = 0;
71  liste->isorder = 0;
72  liste->array = nullptr;
73 
74  CSRList_Realloc(liste, n);
75  return (liste);
76 }
77 
78 static void CSRList_Delete(CSRList_T *liste)
79 {
80  if(liste != nullptr) {
81  free(liste->array);
82  free(liste);
83  }
84 }
85 
86 void CSRList_Add(CSRList_T *liste, const void *data)
87 {
88  liste->n++;
89  CSRList_Realloc(liste, liste->n);
90  liste->isorder = 0;
91  memcpy(&liste->array[(liste->n - 1) * liste->size], data, liste->size);
92 }
93 
94 int CSRList_Nbr(CSRList_T *liste) { return (liste->n); }
95 
96 template <> void linearSystemCSR<double>::preAllocateEntries()
97 {
98  if(_entriesPreAllocated) return;
99  if(_sparsity.getNbRows() == 0) return;
100  INDEX_TYPE nnz = 0;
101  int nbRows = _b->size();
102  for(int i = 0; i < nbRows; i++) {
103  int nInRow;
104  _sparsity.getRow(i, nInRow);
105  nnz += nInRow;
106  }
107  CSRList_Resize_strict(_ai, nnz);
108  CSRList_Resize_strict(_ptr, nnz);
109  INDEX_TYPE *jptr = (INDEX_TYPE *)_jptr->array;
110  INDEX_TYPE *ai = (INDEX_TYPE *)_ai->array;
111  INDEX_TYPE *ptr = (INDEX_TYPE *)_ptr->array;
112  jptr[0] = 0;
113  nnz = 0;
114  for(int i = 0; i < nbRows; i++) {
115  int nInRow;
116  const int *row = _sparsity.getRow(i, nInRow);
117  for(int j = 0; j < nInRow; j++) {
118  ai[nnz] = row[j];
119  ptr[nnz] = nnz + 1;
120  nnz++;
121  }
122  if(nInRow != 0) ptr[nnz - 1] = 0;
123  jptr[i + 1] = nnz;
124  something[i] = (nInRow == 0 ? 0 : 1);
125  }
126  _entriesPreAllocated = true;
127  sorted = true;
128  _sparsity.clear();
129  // we do this after _sparsity.clear so that the peak memory usage is reduced
130  CSRList_Resize_strict(_a, nnz);
131  double *a = (double *)_a->array;
132  for(int i = 0; i < nnz; i++) {
133  a[i] = 0.;
134  }
135 }
136 
137 template <> void linearSystemCSR<std::complex<double> >::preAllocateEntries()
138 {
139  if(_entriesPreAllocated) return;
140  if(_sparsity.getNbRows() == 0) return;
141  INDEX_TYPE nnz = 0;
142  int nbRows = _b->size();
143  for(int i = 0; i < nbRows; i++) {
144  int nInRow;
145  _sparsity.getRow(i, nInRow);
146  nnz += nInRow;
147  }
148  CSRList_Resize_strict(_ai, nnz);
149  CSRList_Resize_strict(_ptr, nnz);
150  INDEX_TYPE *jptr = (INDEX_TYPE *)_jptr->array;
151  INDEX_TYPE *ai = (INDEX_TYPE *)_ai->array;
152  INDEX_TYPE *ptr = (INDEX_TYPE *)_ptr->array;
153  jptr[0] = 0;
154  nnz = 0;
155  for(int i = 0; i < nbRows; i++) {
156  int nInRow;
157  const int *row = _sparsity.getRow(i, nInRow);
158  for(int j = 0; j < nInRow; j++) {
159  ai[nnz] = row[j];
160  ptr[nnz] = nnz + 1;
161  nnz++;
162  }
163  if(nInRow != 0) ptr[nnz - 1] = 0;
164  jptr[i + 1] = nnz;
165  something[i] = (nInRow == 0 ? 0 : 1);
166  }
167  _entriesPreAllocated = true;
168  sorted = true;
169  _sparsity.clear();
170  // we do this after _sparsity.clear so that the peak memory usage is reduced
171  CSRList_Resize_strict(_a, nnz);
172  std::complex<double> *a = (std::complex<double> *)_a->array;
173  for(int i = 0; i < nnz; i++) {
174  a[i] = std::complex<double>();
175  }
176 }
177 
178 template <> void linearSystemCSR<double>::allocate(int nbRows)
179 {
180  if(_a) {
181  CSRList_Delete(_a);
182  CSRList_Delete(_ai);
183  CSRList_Delete(_ptr);
184  CSRList_Delete(_jptr);
185  delete _x;
186  delete _b;
187  delete[] something;
188  }
189 
190  if(nbRows == 0) {
191  _a = nullptr;
192  _ai = nullptr;
193  _ptr = nullptr;
194  _jptr = nullptr;
195  _b = nullptr;
196  _x = nullptr;
197  sorted = false;
198  something = nullptr;
199  return;
200  }
201 
202  _a = CSRList_Create(nbRows, nbRows, sizeof(double));
203  _ai = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
204  _ptr = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
205  _jptr = CSRList_Create(nbRows + 1, nbRows, sizeof(INDEX_TYPE));
206 
207  something = new char[nbRows];
208 
209  for(int i = 0; i < nbRows; i++) something[i] = 0;
210 
211  _b = new std::vector<double>(nbRows);
212  _x = new std::vector<double>(nbRows);
213 }
214 
215 template <> void linearSystemCSR<std::complex<double> >::allocate(int nbRows)
216 {
217  if(_a) {
218  CSRList_Delete(_a);
219  CSRList_Delete(_ai);
220  CSRList_Delete(_ptr);
221  CSRList_Delete(_jptr);
222  delete _x;
223  delete _b;
224  delete[] something;
225  }
226 
227  if(nbRows == 0) {
228  _a = nullptr;
229  _ai = nullptr;
230  _ptr = nullptr;
231  _jptr = nullptr;
232  _b = nullptr;
233  _x = nullptr;
234  sorted = false;
235  something = nullptr;
236  return;
237  }
238 
239  _a = CSRList_Create(nbRows, nbRows, sizeof(std::complex<double>));
240  _ai = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
241  _ptr = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
242  _jptr = CSRList_Create(nbRows + 1, nbRows, sizeof(INDEX_TYPE));
243 
244  something = new char[nbRows];
245 
246  for(int i = 0; i < nbRows; i++) something[i] = 0;
247 
248  _b = new std::vector<std::complex<double> >(nbRows);
249  _x = new std::vector<std::complex<double> >(nbRows);
250 }
251 
252 const int NSTACK = 50;
253 const unsigned int M_sort2 = 7;
254 
255 static void free_ivector(int *v, long nl, long nh)
256 {
257  // free an int vector allocated with ivector()
258  free((char *)(v + nl - 1));
259 }
260 
261 static int *ivector(long nl, long nh)
262 {
263  // allocate an int vector with subscript range v[nl..nh]
264  int *v;
265  v = (int *)malloc((size_t)((nh - nl + 2) * sizeof(int)));
266  if(!v) fprintf(stderr, "allocation failure in ivector()\n");
267  return v - nl + 1;
268 }
269 
270 static int cmpij(INDEX_TYPE ai, INDEX_TYPE aj, INDEX_TYPE bi, INDEX_TYPE bj)
271 {
272  if(ai < bi) return -1;
273  if(ai > bi) return 1;
274  if(aj < bj) return -1;
275  if(aj > bj) return 1;
276  return 0;
277 }
278 
279 template <class scalar>
280 static void _sort2_xkws(unsigned long n, scalar arr[], INDEX_TYPE ai[],
281  INDEX_TYPE aj[])
282 {
283  unsigned long i, ir = n, j, k, l = 1;
284  int *istack, jstack = 0;
285  INDEX_TYPE tempi;
286  scalar a, temp;
287  int b, c;
288 
289  istack = ivector(1, NSTACK);
290  for(;;) {
291  if(ir - l < M_sort2) {
292  for(j = l + 1; j <= ir; j++) {
293  a = arr[j - 1];
294  b = ai[j - 1];
295  c = aj[j - 1];
296  for(i = j - 1; i >= 1; i--) {
297  if(cmpij(ai[i - 1], aj[i - 1], b, c) <= 0) break;
298  arr[i + 1 - 1] = arr[i - 1];
299  ai[i + 1 - 1] = ai[i - 1];
300  aj[i + 1 - 1] = aj[i - 1];
301  }
302  arr[i + 1 - 1] = a;
303  ai[i + 1 - 1] = b;
304  aj[i + 1 - 1] = c;
305  }
306  if(!jstack) {
307  free_ivector(istack, 1, NSTACK);
308  return;
309  }
310  ir = istack[jstack];
311  l = istack[jstack - 1];
312  jstack -= 2;
313  }
314  else {
315  k = (l + ir) >> 1;
316  SWAP(arr[k - 1], arr[l + 1 - 1])
317  SWAPI(ai[k - 1], ai[l + 1 - 1])
318  SWAPI(aj[k - 1], aj[l + 1 - 1])
319  if(cmpij(ai[l + 1 - 1], aj[l + 1 - 1], ai[ir - 1], aj[ir - 1]) > 0) {
320  SWAP(arr[l + 1 - 1], arr[ir - 1])
321  SWAPI(ai[l + 1 - 1], ai[ir - 1])
322  SWAPI(aj[l + 1 - 1], aj[ir - 1])
323  }
324  if(cmpij(ai[l - 1], aj[l - 1], ai[ir - 1], aj[ir - 1]) > 0) {
325  SWAP(arr[l - 1], arr[ir - 1])
326  SWAPI(ai[l - 1], ai[ir - 1])
327  SWAPI(aj[l - 1], aj[ir - 1])
328  }
329  if(cmpij(ai[l + 1 - 1], aj[l + 1 - 1], ai[l - 1], aj[l - 1]) > 0) {
330  SWAP(arr[l + 1 - 1], arr[l - 1])
331  SWAPI(ai[l + 1 - 1], ai[l - 1])
332  SWAPI(aj[l + 1 - 1], aj[l - 1])
333  }
334  i = l + 1;
335  j = ir;
336  a = arr[l - 1];
337  b = ai[l - 1];
338  c = aj[l - 1];
339  for(;;) {
340  do
341  i++;
342  while(cmpij(ai[i - 1], aj[i - 1], b, c) < 0);
343  do
344  j--;
345  while(cmpij(ai[j - 1], aj[j - 1], b, c) > 0);
346  if(j < i) break;
347  SWAP(arr[i - 1], arr[j - 1])
348  SWAPI(ai[i - 1], ai[j - 1])
349  SWAPI(aj[i - 1], aj[j - 1])
350  }
351  arr[l - 1] = arr[j - 1];
352  arr[j - 1] = a;
353  ai[l - 1] = ai[j - 1];
354  ai[j - 1] = b;
355  aj[l - 1] = aj[j - 1];
356  aj[j - 1] = c;
357  jstack += 2;
358  if(jstack > NSTACK) {
359  Msg::Error("NSTACK too small while sorting the columns of the matrix");
360  return;
361  }
362  if(ir - i + 1 >= j - l) {
363  istack[jstack] = ir;
364  istack[jstack - 1] = i;
365  ir = j - 1;
366  }
367  else {
368  istack[jstack] = j - 1;
369  istack[jstack - 1] = l;
370  l = i;
371  }
372  }
373  }
374 }
375 
376 template <class scalar>
377 void sortColumns_(int NbLines, int nnz, INDEX_TYPE *ptr, INDEX_TYPE *jptr,
378  INDEX_TYPE *ai, scalar *a)
379 {
380  // replace pointers by lines
381  int *count = new int[NbLines];
382 
383  for(int i = 0; i < NbLines; i++) {
384  count[i] = 0;
385  INDEX_TYPE _position = jptr[i];
386  while(1) {
387  count[i]++;
388  INDEX_TYPE _position_temp = _position;
389  _position = ptr[_position];
390  ptr[_position_temp] = i;
391  if(_position == 0) break;
392  }
393  }
394  _sort2_xkws<scalar>(nnz, a, ptr, ai);
395  jptr[0] = 0;
396  for(int i = 1; i <= NbLines; i++) {
397  jptr[i] = jptr[i - 1] + count[i - 1];
398  }
399 
400  for(int i = 0; i < NbLines; i++) {
401  for(int j = jptr[i]; j < jptr[i + 1] - 1; j++) {
402  ptr[j] = j + 1;
403  }
404  if(jptr[i + 1] != jptr[i]) ptr[jptr[i + 1] - 1] = 0;
405  }
406 
407  delete[] count;
408 }
409 
410 template <>
411 void linearSystemCSR<double>::getMatrix(INDEX_TYPE *&jptr, INDEX_TYPE *&ai,
412  double *&a)
413 {
414  jptr = (INDEX_TYPE *)_jptr->array;
415  ai = (INDEX_TYPE *)_ai->array;
416  a = (double *)_a->array;
417  if(!sorted)
418  sortColumns_(_b->size(), CSRList_Nbr(_a), (INDEX_TYPE *)_ptr->array, jptr,
419  ai, a);
420  sorted = true;
421 }
422 
423 template <>
424 void linearSystemCSR<std::complex<double> >::getMatrix(INDEX_TYPE *&jptr,
425  INDEX_TYPE *&ai,
426  double *&a)
427 {
428  jptr = (INDEX_TYPE *)_jptr->array;
429  ai = (INDEX_TYPE *)_ai->array;
430  a = (double *)_a->array;
431  if(!sorted)
432  sortColumns_(_b->size(), CSRList_Nbr(_a), (INDEX_TYPE *)_ptr->array, jptr,
433  ai, (std::complex<double> *)a);
434  sorted = true;
435 }
436 
437 #if defined(HAVE_GMM)
438 
439 #undef BB // can be defined by FlGui.h, and clashes with gmm arg name
440 #include <gmm.h>
441 
442 template <> int linearSystemCSRGmm<double>::systemSolve()
443 {
444  if(!sorted)
445  sortColumns_(_b->size(), CSRList_Nbr(_a), (INDEX_TYPE *)_ptr->array,
446  (INDEX_TYPE *)_jptr->array, (INDEX_TYPE *)_ai->array,
447  (double *)_a->array);
448  sorted = true;
449 
450  gmm::csr_matrix_ref<double *, INDEX_TYPE *, INDEX_TYPE *, 0> ref(
451  (double *)_a->array, (INDEX_TYPE *)_ai->array, (INDEX_TYPE *)_jptr->array,
452  _b->size(), _b->size());
453  gmm::csr_matrix<double> M;
454  M.init_with(ref);
455 
456  //gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
457  gmm::ilu_precond<gmm::csr_matrix<double> > P(M);
458  gmm::iteration iter(_tol);
459  iter.set_noisy(_noisy);
460  if(_method == "gmres")
461  gmm::gmres(M, *_x, *_b, P, 100, iter);
462  else
463  gmm::cg(M, *_x, *_b, P, iter);
464  if(!iter.converged())
465  Msg::Warning("Iterative linear solver has not converged (res = %g)",
466  iter.get_res());
467  return 1;
468 }
469 
470 #endif
M_sort2
const unsigned int M_sort2
Definition: linearSystemCSR.cpp:253
CSRList_Delete
static void CSRList_Delete(CSRList_T *liste)
Definition: linearSystemCSR.cpp:78
linearSystemCSR
Definition: linearSystemCSR.h:29
CSRList_T::nmax
int nmax
Definition: linearSystemCSR.h:18
OS.h
CSRList_Create
static CSRList_T * CSRList_Create(int n, int incr, int size)
Definition: linearSystemCSR.cpp:58
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
sortColumns_
void sortColumns_(int NbLines, int nnz, INDEX_TYPE *ptr, INDEX_TYPE *jptr, INDEX_TYPE *ai, scalar *a)
Definition: linearSystemCSR.cpp:377
linearSystemCSR::allocate
virtual void allocate(int)
SWAP
#define SWAP(a, b)
Definition: linearSystemCSR.cpp:15
linearSystemCSR.h
CSRRealloc
static void * CSRRealloc(void *ptr, size_t size)
Definition: linearSystemCSR.cpp:32
GmshMessage.h
CSRList_T::isorder
int isorder
Definition: linearSystemCSR.h:22
nanoflann::allocate
T * allocate(size_t count=1)
Definition: nanoflann.hpp:442
CSRList_T::n
int n
Definition: linearSystemCSR.h:21
CSRList_Add
void CSRList_Add(CSRList_T *liste, const void *data)
Definition: linearSystemCSR.cpp:86
INDEX_TYPE
int INDEX_TYPE
Definition: linearSystemCSR.h:16
CSRMalloc
static void * CSRMalloc(size_t size)
Definition: linearSystemCSR.cpp:24
_sort2_xkws
static void _sort2_xkws(unsigned long n, scalar arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
Definition: linearSystemCSR.cpp:280
CSRList_T::array
char * array
Definition: linearSystemCSR.h:23
linearSystemCSR::getMatrix
virtual void getMatrix(INDEX_TYPE *&jptr, INDEX_TYPE *&ai, double *&a)
CSRList_Nbr
int CSRList_Nbr(CSRList_T *liste)
Definition: linearSystemCSR.cpp:94
ivector
static int * ivector(long nl, long nh)
Definition: linearSystemCSR.cpp:261
CSRList_Realloc
static void CSRList_Realloc(CSRList_T *liste, int n)
Definition: linearSystemCSR.cpp:34
CSRList_T::size
int size
Definition: linearSystemCSR.h:19
CSRList_Resize_strict
static void CSRList_Resize_strict(CSRList_T *liste, int n)
Definition: linearSystemCSR.cpp:51
linearSystemCSR::preAllocateEntries
virtual void preAllocateEntries()
NSTACK
const int NSTACK
Definition: linearSystemCSR.cpp:252
free_ivector
static void free_ivector(int *v, long nl, long nh)
Definition: linearSystemCSR.cpp:255
SWAPI
#define SWAPI(a, b)
Definition: linearSystemCSR.cpp:19
linearSystemCSRGmm::systemSolve
virtual int systemSolve()
Definition: linearSystemCSR.h:190
cmpij
static int cmpij(INDEX_TYPE ai, INDEX_TYPE aj, INDEX_TYPE bi, INDEX_TYPE bj)
Definition: linearSystemCSR.cpp:270
CSRList_T
Definition: linearSystemCSR.h:17
CSRList_T::incr
int incr
Definition: linearSystemCSR.h:20