Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
ordering.c
Go to the documentation of this file.
1 
6 #include "fasp.h"
7 
8 static void dSwapping(REAL *w, const INT i, const INT j);
9 static void iSwapping(INT *w, const INT i, const INT j);
10 static void CMK_ordering (const dCSRmat *, INT, INT, INT, INT, INT *, INT *);
11 
12 /*---------------------------------*/
13 /*-- Public Functions --*/
14 /*---------------------------------*/
15 
31  const INT value,
32  const INT nlist)
33 {
34  INT not_found = 1;
35  INT low, high, m;
36 
37  low = 0;
38  high = nlist - 1;
39  while (not_found && low <= high)
40  {
41  m = (low + high) / 2;
42  if (value < list[m])
43  {
44  high = m - 1;
45  }
46  else if (value > list[m])
47  {
48  low = m + 1;
49  }
50  else
51  {
52  not_found = 0;
53  return m;
54  }
55  }
56 
57  return -1;
58 }
59 
75 INT fasp_aux_unique (INT numbers[],
76  const INT size)
77 {
78  INT i, newsize;
79 
80  if (size==0) return(0);
81 
82  for (newsize=0, i=1; i<size; ++i) {
83  if (numbers[newsize] < numbers[i]) {
84  newsize++;
85  numbers[newsize] = numbers[i];
86  }
87  }
88 
89  return(newsize+1);
90 }
91 
108 void fasp_aux_merge (INT numbers[],
109  INT work[],
110  INT left,
111  INT mid,
112  INT right)
113 {
114  INT i, left_end, num_elements, tmp_pos;
115 
116  left_end = mid - 1;
117  tmp_pos = left;
118  num_elements = right - left + 1;
119 
120  while ((left <= left_end) && (mid <= right)) {
121 
122  if (numbers[left] <= numbers[mid]) // first branch <=
123  {
124  work[tmp_pos] = numbers[left];
125  tmp_pos = tmp_pos + 1;
126  left = left +1;
127  }
128  else // second branch >
129  {
130  work[tmp_pos] = numbers[mid];
131  tmp_pos = tmp_pos + 1;
132  mid = mid + 1;
133  }
134  }
135 
136  while (left <= left_end) {
137  work[tmp_pos] = numbers[left];
138  left = left + 1;
139  tmp_pos = tmp_pos + 1;
140  }
141 
142  while (mid <= right) {
143  work[tmp_pos] = numbers[mid];
144  mid = mid + 1;
145  tmp_pos = tmp_pos + 1;
146  }
147 
148  for (i = 0; i < num_elements; ++i) {
149  numbers[right] = work[right];
150  right = right - 1;
151  }
152 
153 }
154 
170 void fasp_aux_msort (INT numbers[],
171  INT work[],
172  INT left,
173  INT right)
174 {
175  INT mid;
176 
177  if (right > left) {
178  mid = (right + left) / 2;
179  fasp_aux_msort(numbers, work, left, mid);
180  fasp_aux_msort(numbers, work, mid+1, right);
181  fasp_aux_merge(numbers, work, left, mid+1, right);
182  }
183 
184 }
185 
202  INT left,
203  INT right)
204 {
205  INT i, last;
206 
207  if (left >= right) return;
208 
209  iSwapping(a, left, (left+right)/2);
210 
211  last = left;
212  for (i = left+1; i <= right; ++i) {
213  if (a[i] < a[left]) {
214  iSwapping(a, ++last, i);
215  }
216  }
217 
218  iSwapping(a, left, last);
219 
220  fasp_aux_iQuickSort(a, left, last-1);
221  fasp_aux_iQuickSort(a, last+1, right);
222 }
223 
240  INT left,
241  INT right)
242 {
243  INT i, last;
244 
245  if (left >= right) return;
246 
247  dSwapping(a, left, (left+right)/2);
248 
249  last = left;
250  for (i = left+1; i <= right; ++i) {
251  if (a[i] < a[left]) {
252  dSwapping(a, ++last, i);
253  }
254  }
255 
256  dSwapping(a, left, last);
257 
258  fasp_aux_dQuickSort(a, left, last-1);
259  fasp_aux_dQuickSort(a, last+1, right);
260 }
261 
280  INT left,
281  INT right,
282  INT *index)
283 {
284  INT i, last;
285 
286  if (left >= right) return;
287 
288  iSwapping(index, left, (left+right)/2);
289 
290  last = left;
291  for (i = left+1; i <= right; ++i) {
292  if (a[index[i]] < a[index[left]]) {
293  iSwapping(index, ++last, i);
294  }
295  }
296 
297  iSwapping(index, left, last);
298 
299  fasp_aux_iQuickSortIndex(a, left, last-1, index);
300  fasp_aux_iQuickSortIndex(a, last+1, right, index);
301 }
302 
321  INT left,
322  INT right,
323  INT *index)
324 {
325  INT i, last;
326 
327  if (left >= right) return;
328 
329  iSwapping(index, left, (left+right)/2);
330 
331  last = left;
332  for (i = left+1; i <= right; ++i) {
333  if (a[index[i]] < a[index[left]]) {
334  iSwapping(index, ++last, i);
335  }
336  }
337 
338  iSwapping(index, left, last);
339 
340  fasp_aux_dQuickSortIndex(a, left, last-1, index);
341  fasp_aux_dQuickSortIndex(a, last+1, right, index);
342 }
343 
357  INT *order,
358  INT *oindex)
359 {
360  const INT *ia = A->IA;
361  const INT row= A->row;
362 
363  INT i, loc, s, vt, mindg, innz;
364 
365  s = 0;
366  vt = 0;
367  mindg = row+1;
368 
369  // select node with minimal degree
370  for (i=0; i<row; ++i) {
371  innz = ia[i+1] - ia[i];
372  if (innz > 1) {
373  oindex[i] = -innz;
374  if (innz < mindg) {
375  mindg = innz;
376  vt = i;
377  }
378  }
379  else { // order those diagonal rows first
380  oindex[i] = s;
381  order[s] = i;
382  s ++;
383  }
384  }
385 
386  loc = s;
387 
388  // start to order
389  CMK_ordering (A, loc, s, vt, mindg, oindex, order);
390 }
391 
406  INT *order,
407  INT *oindex,
408  INT *rorder)
409 {
410  INT i;
411  INT row = A->row;
412 
413  // Form CMK order
414  fasp_dcsr_CMK_order(A, order, oindex);
415 
416  // Reverse CMK order
417  for (i=0; i<row; ++i) rorder[i] = order[row-1-i];
418 }
419 
420 /*---------------------------------*/
421 /*-- Private Functions --*/
422 /*---------------------------------*/
423 
436 static void iSwapping (INT *w,
437  const INT i,
438  const INT j)
439 {
440  INT temp = w[i];
441  w[i] = w[j];
442  w[j] = temp;
443 }
444 
457 static void dSwapping (REAL *w,
458  const INT i,
459  const INT j)
460 {
461  REAL temp = w[i];
462  w[i] = w[j];
463  w[j] = temp;
464 }
465 
483 static void CMK_ordering (const dCSRmat *A,
484  INT loc,
485  INT s,
486  INT jj,
487  INT mindg,
488  INT *oindex,
489  INT *order)
490 {
491  const INT *ia = A->IA;
492  const INT *ja = A->JA;
493  const INT row= A->row;
494 
495  INT i, j, sp1, k, flag;
496 
497  if (s < row) {
498  order[s] = jj;
499  oindex[jj] = s;
500  }
501 
502  while (loc <= s && s < row) {
503  i = order[loc];
504  sp1 = s+1;
505  // neighbor nodes are priority.
506  for (j=ia[i]+1; j<ia[i+1]; ++j) {
507  k = ja[j];
508  if (oindex[k] < 0){
509  s++;
510  order[s] = k;
511  }
512  }
513  // ordering neighbor nodes by increasing degree
514  if (s > sp1) {
515  while (flag) {
516  flag = 0;
517  for (i=sp1+1; i<=s; ++i) {
518  if (oindex[order[i]] > oindex[order[i-1]]) {
519  j = order[i];
520  order[i] = order[i-1];
521  order[i-1] = j;
522  flag = 1;
523  }
524  }
525  }
526  }
527 
528  for (i=sp1; i<=s; ++i) oindex[order[i]] = i;
529 
530  loc ++;
531  }
532 
533  // deal with remainder
534  if (s < row) {
535  jj = 0;
536  i = 0;
537  while (jj == 0) {
538  i ++;
539  if (i >= row) {
540  mindg++;
541  i = 0;
542  }
543  if (oindex[i] < 0 && (ia[i+1]-ia[i] == mindg)) {
544  jj = i;
545  }
546  }
547 
548  s ++;
549 
550  CMK_ordering (A, loc, s, jj, mindg, oindex, order);
551  }
552 }
553 
554 /*---------------------------------*/
555 /*-- End of File --*/
556 /*---------------------------------*/
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
void fasp_aux_iQuickSortIndex(INT *a, INT left, INT right, INT *index)
Reorder the index of (INT type) so that 'a' is in ascending order.
Definition: ordering.c:279
void fasp_aux_msort(INT numbers[], INT work[], INT left, INT right)
Sort the INT array in ascending order with the merge sort algorithm.
Definition: ordering.c:170
void fasp_dcsr_CMK_order(const dCSRmat *A, INT *order, INT *oindex)
Ordering vertices of matrix graph corresponding to A.
Definition: ordering.c:356
void fasp_aux_dQuickSortIndex(REAL *a, INT left, INT right, INT *index)
Reorder the index of (REAL type) so that 'a' is ascending in such order.
Definition: ordering.c:320
#define INT
Definition: fasp.h:64
INT fasp_BinarySearch(INT *list, const INT value, const INT nlist)
Binary Search.
Definition: ordering.c:30
Main header file for FASP.
void fasp_aux_merge(INT numbers[], INT work[], INT left, INT mid, INT right)
Merge two sorted arrays.
Definition: ordering.c:108
void fasp_dcsr_RCMK_order(const dCSRmat *A, INT *order, INT *oindex, INT *rorder)
Resverse CMK ordering.
Definition: ordering.c:405
void fasp_aux_dQuickSort(REAL *a, INT left, INT right)
Sort the array (REAL type) in ascending order with the quick sorting algorithm.
Definition: ordering.c:239
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT fasp_aux_unique(INT numbers[], const INT size)
Remove duplicates in an sorted (ascending order) array.
Definition: ordering.c:75
void fasp_aux_iQuickSort(INT *a, INT left, INT right)
Sort the array (INT type) in ascending order with the quick sorting algorithm.
Definition: ordering.c:201
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160