Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
coarsening_cr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 
15 static INT GraphAdd(Link *list, INT *head, INT *tail, INT index, INT istack);
16 static INT GraphRemove(Link *list, INT *head, INT *tail, INT index);
17 static INT indset(INT cand, INT cpt, INT fpt, INT *ia, INT *ja, INT n, INT *cf, REAL *ma);
18 
19 /*---------------------------------*/
20 /*-- Public Functions --*/
21 /*---------------------------------*/
22 
43  const INT i_n,
44  dCSRmat *A,
45  ivector *vertices,
46  AMG_param *param)
47 {
48  const SHORT prtlvl = param->print_level;
49 
50  // local variables
51  INT cand=0,cpt=-1,fpt=1; // internal labeling
52  INT nc,ns=1; // # cpts, # stages
53  INT i,j,in1,nu=3,num1 = nu-1; // nu is number of cr sweeps
54  INT *cf=NULL,*ia=NULL,*ja=NULL;
55 
56  REAL temp0=0.0e0,temp1=0.0e0,rho=0.0e0,tg=8.0e-01;
57  REAL *a=NULL;
58 
59  /* WORKING MEMORY -- b not needed, remove later */
60  REAL *b=NULL,*u=NULL,*ma=NULL;
61 
62 #ifdef _OPENMP
63  // variables for OpenMP
64  INT myid, mybegin, myend;
65  REAL sub_temp0 = 0.;
66  INT nthreads = FASP_GET_NUM_THREADS();
67 #endif
68 
69  ia = A->IA;
70  ja = A->JA;
71  a = A->val;
72  cf = vertices->val;
73  if (i_0 == 0) {
74  in1 = i_n+1;
75  } else {
76  in1 = i_n;
77  }
78 
79  /* CF, RHS, INITIAL GUESS, and MEAS. ARRAY */
80  cf = (INT*)fasp_mem_calloc(in1,sizeof(INT));
81  b = (REAL*)fasp_mem_calloc(in1,sizeof(REAL));
82  u = (REAL*)fasp_mem_calloc(in1,sizeof(REAL));
83  ma = (REAL*)fasp_mem_calloc(in1,sizeof(REAL));
84 
85 #ifdef _OPENMP
86 #pragma omp parallel for if(i_n>OPENMP_HOLDS)
87 #endif
88  for(i=i_0;i<=i_n;++i) {
89  b[i] = 0.0e0; // ZERO RHS
90  cf[i] = fpt; // ALL FPTS
91  }
92 
94  while(1) {
95  nc = 0;
96 #ifdef _OPENMP
97 #pragma omp parallel for if(i_n>OPENMP_HOLDS)
98 #endif
99  for(i=i_0; i<=i_n; ++i) {
100  if (cf[i] == cpt) {
101  nc += 1;
102  u[i] = 0.0e0;
103  } else {
104  u[i] = 1.0e0;
105  }
106  }
107 #ifdef _OPENMP
108 #pragma omp parallel for private(myid,mybegin,myend,i,j,sub_temp0) if(nu>OPENMP_HOLDS)
109  for (myid=0; myid<nthreads; ++myid) {
110  FASP_GET_START_END(myid, nthreads, nu-i_0+1, &mybegin, &myend);
111  mybegin += i_0; myend += i_0;
112  for (i=mybegin; i<myend; ++i) {
113 #else
114  for (i=i_0;i<=nu;++i) {
115 #endif
116  if (i == num1)
117  for (j = i_0; j<= i_n; ++j) {
118  if (cf[j] == fpt) {
119 #ifdef _OPENMP
120  sub_temp0 += u[j]*u[j];
121 #else
122  temp0 += u[j]*u[j];
123 #endif
124  }
125  }
126  fasp_smoother_dcsr_gscr(fpt,i_n,u,ia,ja,a,b,1,cf);
127  }
128 #ifdef _OPENMP
129 #pragma omp critical(temp0)
130  temp0 += sub_temp0;
131  }
132 #endif
133 
134 #ifdef _OPENMP
135 #pragma omp parallel for reduction(+:temp1) if(i_n>OPENMP_HOLDS)
136 #endif
137  for (i = i_0; i<= i_n; ++i) {
138  if (cf[i] == fpt) {
139  temp1 += u[i]*u[i];
140  }
141  }
142  rho = sqrt(temp1)/sqrt(temp0);
143  temp0 = 0.0e0;
144  temp1 = 0.0e0;
145 
146  if ( prtlvl > PRINT_MIN ) printf("rho=%2.13lf\n",rho);
147 
148  if ( rho > tg ) {
149  /* FORM CAND. SET & COMPUTE IND SET */
150  temp0 = 0.0e0;
151 
152  for (i = i_0; i<= i_n; ++i) {
153  temp1 = fabs(u[i]);
154  if (cf[i] == cpt && temp1 > 0.0e0) {
155  temp0 = temp1; // max.
156  }
157  }
158  temp1 = 0.0e0;
159  if (ns == 1) {
160  temp1 = pow(0.3, nu);
161  } else {
162  temp1 = 0.5;
163  }
164 
165 #ifdef _OPENMP
166 #pragma omp parallel for if(i_n>OPENMP_HOLDS)
167 #endif
168  for (i = i_0; i <= i_n; ++i) {
169  if (cf[i] == fpt && fabs(u[i])/temp0 > temp1 && ia[i+1]-ia[i] > 1)
170  cf[i] = cand;
171  }
172  temp1 = 0.0e0;
173  indset(cand,cpt,fpt,ia,ja,i_n,cf,ma);
174  ns++;
175  } else {
176  /* back to fasp labeling */
177 
178 #ifdef _OPENMP
179 #pragma omp parallel for if(i_n>OPENMP_HOLDS)
180 #endif
181  for (i = i_0; i<= i_n; ++i) {
182  if (cf[i] == cpt) {
183  cf[i] = 1; // cpt
184  } else {
185  cf[i] = 0; // fpt
186  }
187  // printf("cf[%i] = %i\n",i,cf[i]);
188  }
189  vertices->row=i_n;
190  if ( prtlvl >= PRINT_MORE ) printf("vertices = %i\n",vertices->row);
191  vertices->val= cf;
192  if ( prtlvl >= PRINT_MORE ) printf("nc=%i\n",nc);
193  break;
194  }
195  }
196 
197  fasp_mem_free(u);
198  fasp_mem_free(b);
199  fasp_mem_free(ma);
200 
201  return nc;
202 }
203 
204 /*---------------------------------*/
205 /*-- Private Functions --*/
206 /*---------------------------------*/
207 
219 static INT GraphAdd(Link *list, INT *head, INT *tail, INT index, INT istack)
220 {
221  INT prev = tail[-istack];
222 
223  list[index].prev = prev;
224  if (prev < 0)
225  head[-istack] = index;
226  else
227  list[prev].next = index;
228  list[index].next = -istack;
229  tail[-istack] = index;
230 
231  return 0;
232 }
233 
244 static INT GraphRemove(Link *list, INT *head, INT *tail, INT index)
245 {
246  INT prev = list[index].prev;
247  INT next = list[index].next;
248 
249  if (prev < 0)
250  head[prev] = next;
251  else
252  list[prev].next = next;
253  if (next < 0)
254  tail[next] = prev;
255  else
256  list[next].prev = prev;
257 
258  return 0;
259 }
260 
280 static INT indset(INT cand, INT cpt, INT fpt, INT *ia, INT *ja, INT n, INT *cf, REAL *ma)
281 {
282  /* ma: candidates >= 1, cpts = -1, otherwise = 0
283  * Note: graph contains candidates only */
284 
285  Link *list;
286  INT *head, *head_mem;
287  INT *tail, *tail_mem;
288 
289  INT i, ji, jj, jl, index, istack, stack_size;
290 
291 #ifdef _OPENMP
292  // variables for OpenMP
293  INT myid, mybegin, myend;
294  INT sub_istack = 0;
295  INT nthreads = FASP_GET_NUM_THREADS();
296 #endif
297 
298  istack = 0;
299 
300 #ifdef _OPENMP
301 #pragma omp parallel for private(myid,mybegin,myend,i,ji,jj,sub_istack) if(n>OPENMP_HOLDS)
302  for (myid=0; myid<nthreads; ++myid) {
303  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
304  for (i=mybegin; i<myend; ++i) {
305 #else
306  for (i = 0; i < n; ++i) {
307 #endif
308  if (cf[i] == cand) {
309  ma[i] = 1;
310  for (ji = ia[i]+1; ji < ia[i+1]; ++ji) {
311  jj = ja[ji];
312  if (cf[jj] != cpt) {
313  ma[i]++;
314  }
315  }
316 #ifdef _OPENMP
317  if (ma[i] > sub_istack) {
318  sub_istack = (INT) ma[i];
319  }
320 #else
321  if (ma[i] > istack) {
322  istack = (INT) ma[i];
323  }
324 #endif
325  } else if (cf[i] == cpt) {
326  ma[i] = -1;
327  } else {
328  ma[i] = 0;
329  }
330  }
331 #ifdef _OPENMP
332 #pragma omp critical(istack)
333  if (sub_istack > istack) istack = sub_istack;
334  }
335 #endif
336  stack_size = 2*istack;
337 
338  /* INITIALIZE GRAPH */
339  list = (Link*)fasp_mem_calloc(n,sizeof(Link));
340  head_mem = (INT*)fasp_mem_calloc(stack_size,sizeof(INT));
341  tail_mem = (INT*)fasp_mem_calloc(stack_size,sizeof(INT));
342  head = head_mem + stack_size;
343  tail = tail_mem + stack_size;
344 
345 #ifdef _OPENMP
346 #pragma omp parallel for if(stack_size>OPENMP_HOLDS)
347 #endif
348  for (i = -1; i >= -stack_size; i--) {
349  head[i] = i;
350  tail[i] = i;
351  }
352 
353 #ifdef _OPENMP
354 #pragma omp parallel for if(stack_size>OPENMP_HOLDS)
355 #endif
356  for (i = 0; i < n; ++i) {
357  if (ma[i] > 0) {
358  GraphAdd(list, head, tail, i, (INT) ma[i]);
359  }
360  }
361 
362  while (istack > 0) {
363  /* i with maximal measure is at the head of the stacks */
364  i = head[-istack];
365  /* make i a c-point */
366  cf[i] = cpt;
367  ma[i] = -1;
368  /* remove i from graph */
369  GraphRemove(list, head, tail, i);
370 
371  /* update neighbors and neighbors-of-neighbors */
372  for (ji = ia[i]+1; ji < ia[i+1]; ++ji) {
373  jj = ja[ji];
374  /* if not "decided" c or f */
375  if (ma[jj] > -1) {
376  /* if a candidate, remove jj from graph */
377  if (ma[jj] > 0) {
378  GraphRemove(list, head, tail, jj);
379  }
380  /* make jj an f-point and mark "decided" */
381  cf[jj] = fpt;
382  ma[jj] = -1;
383 
384  for (jl = ia[jj]+1; jl < ia[jj+1]; jl++) {
385  index = ja[jl];
386  /* if a candidate, increase likehood of being chosen */
387  if (ma[index] > 0) {
388  ma[index]++;
389  /* move index in graph */
390  GraphRemove(list, head, tail, index);
391  GraphAdd(list, head, tail, index, (INT) ma[index]);
392  if (ma[index] > istack) {
393  istack = (INT) ma[index];
394  }
395  }
396  }
397  }
398  }
399 
400  /* reset istack to point to the biggest non-empty stack */
401  for ( ; istack > 0; istack--) {
402  /* if non-negative, break */
403  if (head[-istack] > -1) {
404  break;
405  }
406  }
407  }
408 
409  fasp_mem_free(list);
410  fasp_mem_free(head_mem);
411  fasp_mem_free(tail_mem);
412 
413  return 0;
414 }
415 
416 /*---------------------------------*/
417 /*-- End of File --*/
418 /*---------------------------------*/
void fasp_smoother_dcsr_gscr(INT pt, INT n, REAL *u, INT *ia, INT *ja, REAL *a, REAL *b, INT L, INT *CF)
Gauss Seidel method restriced to a block.
Parameters for AMG solver.
Definition: fasp.h:583
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
REAL * val
nonzero entries of A
Definition: fasp.h:166
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
#define INT
Definition: fasp.h:64
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
Main header file for FASP.
#define PRINT_MORE
Definition: fasp_const.h:82
INT * val
actual vector entries
Definition: fasp.h:362
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
SHORT print_level
print level for AMG
Definition: fasp.h:589
INT row
number of rows
Definition: fasp.h:359
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
INT fasp_amg_coarsening_cr(const INT i_0, const INT i_n, dCSRmat *A, ivector *vertices, AMG_param *param)
CR coarsening.
Definition: coarsening_cr.c:42
#define PRINT_MIN
Definition: fasp_const.h:80