/[escript]/trunk/paso/src/AMG.c
ViewVC logotype

Contents of /trunk/paso/src/AMG.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3440 - (show annotations)
Fri Jan 14 00:04:53 2011 UTC (8 years, 10 months ago) by gross
File MIME type: text/plain
File size: 22760 byte(s)
a second AMG interpolation implemented
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: AMG preconditioner */
18
19 /**************************************************************/
20
21 /* Author: artak@uq.edu.au, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #define SHOW_TIMING FALSE
26
27 #include "Paso.h"
28 #include "Preconditioner.h"
29 #include "Options.h"
30 #include "PasoUtil.h"
31 #include "UMFPACK.h"
32 #include "MKL.h"
33
34 /**************************************************************/
35
36 /* free all memory used by AMG */
37
38 void Paso_Preconditioner_LocalAMG_free(Paso_Preconditioner_LocalAMG * in) {
39 if (in!=NULL) {
40 Paso_Preconditioner_LocalSmoother_free(in->Smoother);
41 Paso_SparseMatrix_free(in->P);
42 Paso_SparseMatrix_free(in->R);
43 Paso_SparseMatrix_free(in->A_C);
44 Paso_Preconditioner_LocalAMG_free(in->AMG_C);
45 MEMFREE(in->r);
46 MEMFREE(in->x_C);
47 MEMFREE(in->b_C);
48
49
50 MEMFREE(in);
51 }
52 }
53
54 index_t Paso_Preconditioner_LocalAMG_getMaxLevel(const Paso_Preconditioner_LocalAMG * in) {
55 if (in->AMG_C == NULL) {
56 return in->level;
57 } else {
58 return Paso_Preconditioner_LocalAMG_getMaxLevel(in->AMG_C);
59 }
60 }
61 double Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(const Paso_Preconditioner_LocalAMG * in) {
62 if (in->AMG_C == NULL) {
63 if (in->A_C == NULL) {
64 return 1.;
65 } else {
66 return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);
67 }
68 } else {
69 return Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(in->AMG_C);
70 }
71 }
72 dim_t Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(const Paso_Preconditioner_LocalAMG * in) {
73 if (in->AMG_C == NULL) {
74 if (in->A_C == NULL) {
75 return 0;
76 } else {
77 return in->A_C->numRows;
78 }
79 } else {
80 return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);
81 }
82 }
83 /*****************************************************************
84
85 constructs AMG
86
87 ******************************************************************/
88 Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
89
90 Paso_Preconditioner_LocalAMG* out=NULL;
91 bool_t verbose=options->verbose;
92
93 Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;
94 const dim_t n=A_p->numRows;
95 const dim_t n_block=A_p->row_block_size;
96 index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree_S=NULL;
97 dim_t n_F=0, n_C=0, i;
98 double time0=0;
99 const double theta = options->coarsening_threshold;
100 const double tau = options->diagonal_dominance_threshold;
101
102
103 /*
104 is the input matrix A suitable for coarsening
105
106 */
107 if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {
108
109 if (verbose) {
110 /*
111 print stopping condition:
112 - 'SPAR' = min_coarse_matrix_sparsity exceeded
113 - 'SIZE' = min_coarse_matrix_size exceeded
114 - 'LEVEL' = level_max exceeded
115 */
116 printf("Paso_Preconditioner AMG: termination of coarsening by ");
117
118 if (A_p->pattern->len >= options->min_coarse_sparsity * n * n)
119 printf("SPAR ");
120
121 if (n <= options->min_coarse_matrix_size)
122 printf("SIZE ");
123
124 if (level > options->level_max)
125 printf("LEVEL ");
126
127 printf("\n");
128
129 printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
130 level, options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size);
131
132 }
133
134 return NULL;
135 }
136 /* Start Coarsening : */
137
138 F_marker=TMPMEMALLOC(n,index_t);
139 counter=TMPMEMALLOC(n,index_t);
140 degree_S=TMPMEMALLOC(n, dim_t);
141 S=TMPMEMALLOC(A_p->pattern->len, index_t);
142 if ( !( Esys_checkPtr(F_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree_S) || Esys_checkPtr(S) ) ) {
143 /*
144 set splitting of unknows:
145
146 */
147 time0=Esys_timer();
148 if (n_block>1) {
149 Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree_S, S, theta,tau);
150 } else {
151 Paso_Preconditioner_AMG_setStrongConnections(A_p, degree_S, S, theta,tau);
152 }
153 Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree_S, S, F_marker, options->usePanel);
154
155 /* in BoomerAMG interpolation is used FF connectiovity is required :*/
156 if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
157 Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);
158 options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
159
160 if (Esys_noError() ) {
161 #pragma omp parallel for private(i) schedule(static)
162 for (i = 0; i < n; ++i) F_marker[i]=(F_marker[i] == PASO_AMG_IN_F);
163
164 /*
165 count number of unkowns to be eliminated:
166 */
167 n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);
168 n_C=n-n_F;
169 if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
170
171 if ( n_F == 0 ) { /* is a nasty case. a direct solver should be used, return NULL */
172 out = NULL;
173 } else {
174 out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);
175 if (! Esys_checkPtr(out)) {
176 out->level = level;
177 out->n = n;
178 out->n_F = n_F;
179 out->n_block = n_block;
180 out->A_C = NULL;
181 out->P = NULL;
182 out->R = NULL;
183 out->post_sweeps = options->post_sweeps;
184 out->pre_sweeps = options->pre_sweeps;
185 out->r = NULL;
186 out->x_C = NULL;
187 out->b_C = NULL;
188 out->AMG_C = NULL;
189 out->Smoother=NULL;
190 }
191 mask_C=TMPMEMALLOC(n,index_t);
192 rows_in_F=TMPMEMALLOC(n_F,index_t);
193 Esys_checkPtr(mask_C);
194 Esys_checkPtr(rows_in_F);
195 if ( Esys_noError() ) {
196
197 out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
198
199 if (n_C != 0) {
200 /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
201
202 /* allocate helpers :*/
203 out->x_C=MEMALLOC(n_block*n_C,double);
204 out->b_C=MEMALLOC(n_block*n_C,double);
205 out->r=MEMALLOC(n_block*n,double);
206
207 Esys_checkPtr(out->r);
208 Esys_checkPtr(out->x_C);
209 Esys_checkPtr(out->b_C);
210
211 if ( Esys_noError() ) {
212 /* creates index for F:*/
213 #pragma omp parallel private(i)
214 {
215 #pragma omp for schedule(static)
216 for (i = 0; i < n; ++i) {
217 if (F_marker[i]) rows_in_F[counter[i]]=i;
218 }
219 }
220 /* create mask of C nodes with value >-1 gives new id */
221 i=Paso_Util_cumsum_maskedFalse(n,counter, F_marker);
222
223 #pragma omp parallel for private(i) schedule(static)
224 for (i = 0; i < n; ++i) {
225 if (F_marker[i]) {
226 mask_C[i]=-1;
227 } else {
228 mask_C[i]=counter[i];;
229 }
230 }
231 /*
232 get Prolongation :
233 */
234 time0=Esys_timer();
235 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,A_p->pattern->ptr, degree_S,S,n_C,mask_C, options->interpolation_method);
236 if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);
237 }
238 /*
239 construct Restriction operator as transposed of Prolongation operator:
240 */
241 if ( Esys_noError()) {
242 time0=Esys_timer();
243 out->R=Paso_SparseMatrix_getTranspose(out->P);
244 if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
245 }
246 /*
247 construct coarse level matrix:
248 */
249 if ( Esys_noError()) {
250 time0=Esys_timer();
251 Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
252 A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
253 Paso_SparseMatrix_free(Atemp);
254 if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);
255 }
256
257
258 /*
259 constructe courser level:
260
261 */
262 if ( Esys_noError()) {
263 out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);
264 }
265 if ( Esys_noError()) {
266 if ( out->AMG_C == NULL ) {
267 out->reordering = options->reordering;
268 out->refinements = options->coarse_matrix_refinements;
269 /* no coarse level matrix has been constructed. use direct solver */
270 #ifdef MKL
271 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
272 Paso_SparseMatrix_free(A_C);
273 out->A_C->solver_package = PASO_MKL;
274 if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
275 #else
276 #ifdef UMFPACK
277 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
278 Paso_SparseMatrix_free(A_C);
279 out->A_C->solver_package = PASO_UMFPACK;
280 if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
281 #else
282 out->A_C=A_C;
283 out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);
284 out->A_C->solver_package = PASO_SMOOTHER;
285 if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
286 #endif
287 #endif
288 } else {
289 /* finally we set some helpers for the solver step */
290 out->A_C=A_C;
291 }
292 }
293 }
294 }
295 TMPMEMFREE(mask_C);
296 TMPMEMFREE(rows_in_F);
297 }
298 }
299 }
300 TMPMEMFREE(counter);
301 TMPMEMFREE(F_marker);
302 TMPMEMFREE(degree_S);
303 TMPMEMFREE(S);
304
305 if (Esys_noError()) {
306 return out;
307 } else {
308 Paso_Preconditioner_LocalAMG_free(out);
309 return NULL;
310 }
311 }
312
313
314 void Paso_Preconditioner_LocalAMG_solve(Paso_SparseMatrix* A, Paso_Preconditioner_LocalAMG * amg, double * x, double * b) {
315 const dim_t n = amg->n * amg->n_block;
316 double time0=0;
317 const dim_t post_sweeps=amg->post_sweeps;
318 const dim_t pre_sweeps=amg->pre_sweeps;
319
320 /* presmoothing */
321 time0=Esys_timer();
322 Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
323 time0=Esys_timer()-time0;
324 if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);
325 /* end of presmoothing */
326
327 if (amg->n_F < amg->n) { /* is there work on the coarse level? */
328 time0=Esys_timer();
329 Paso_Copy(n, amg->r, b); /* r <- b */
330 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
331 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C); /* b_c = R*r */
332 time0=Esys_timer()-time0;
333 /* coarse level solve */
334 if ( amg->AMG_C == NULL) {
335 time0=Esys_timer();
336 /* A_C is the coarsest level */
337 switch (amg->A_C->solver_package) {
338 case (PASO_MKL):
339 Paso_MKL(amg->A_C, amg->x_C,amg->b_C, amg->reordering, amg->refinements, SHOW_TIMING);
340 break;
341 case (PASO_UMFPACK):
342 Paso_UMFPACK(amg->A_C, amg->x_C,amg->b_C, amg->refinements, SHOW_TIMING);
343 break;
344 case (PASO_SMOOTHER):
345 Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);
346 break;
347 }
348 if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
349 } else {
350 Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C) */
351 }
352 time0=time0+Esys_timer();
353 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */
354
355 /*postsmoothing*/
356
357 /*solve Ax=b with initial guess x */
358 time0=Esys_timer();
359 Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
360 time0=Esys_timer()-time0;
361 if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
362 /*end of postsmoothing*/
363
364 }
365 return;
366 }
367
368 /* theta = threshold for strong connections */
369 /* tau = threshold for diagonal dominance */
370
371 /*S_i={j \in N_i; i strongly coupled to j}
372
373 in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
374 */
375
376 void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,
377 dim_t *degree_S, index_t *S,
378 const double theta, const double tau)
379 {
380 const dim_t n=A->numRows;
381 index_t iptr, i,j;
382 dim_t kdeg;
383 double max_offdiagonal, threshold, sum_row, main_row, fnorm;
384
385
386 #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)
387 for (i=0;i<n;++i) {
388 max_offdiagonal = 0.;
389 sum_row=0;
390 main_row=0;
391 #pragma ivdep
392 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
393 j=A->pattern->index[iptr];
394 fnorm=ABS(A->val[iptr]);
395
396 if( j != i) {
397 max_offdiagonal = MAX(max_offdiagonal,fnorm);
398 sum_row+=fnorm;
399 } else {
400 main_row=fnorm;
401 }
402 }
403 threshold = theta*max_offdiagonal;
404 kdeg=0;
405
406 if (tau*main_row < sum_row) { /* no diagonal domainance */
407 #pragma ivdep
408 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
409 j=A->pattern->index[iptr];
410 if(ABS(A->val[iptr])>threshold && i!=j) {
411 S[A->pattern->ptr[i]+kdeg] = j;
412 kdeg++;
413 }
414 }
415 }
416 degree_S[i]=kdeg;
417 }
418
419 }
420
421 /* theta = threshold for strong connections */
422 /* tau = threshold for diagonal dominance */
423 /*S_i={j \in N_i; i strongly coupled to j}
424
425 in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
426 */
427 void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,
428 dim_t *degree_S, index_t *S,
429 const double theta, const double tau)
430
431 {
432 const dim_t n_block=A->row_block_size;
433 const dim_t n=A->numRows;
434 index_t iptr, i,j, bi;
435 dim_t kdeg, max_deg;
436 register double max_offdiagonal, threshold, fnorm, sum_row, main_row;
437 double *rtmp;
438
439
440 #pragma omp parallel private(i,iptr,max_offdiagonal, kdeg, threshold,j, max_deg, fnorm, sum_row, main_row, rtmp)
441 {
442 max_deg=0;
443 #pragma omp for schedule(static)
444 for (i=0;i<n;++i) max_deg=MAX(max_deg, A->pattern->ptr[i+1]-A->pattern->ptr[i]);
445
446 rtmp=TMPMEMALLOC(max_deg, double);
447
448 #pragma omp for schedule(static)
449 for (i=0;i<n;++i) {
450
451 max_offdiagonal = 0.;
452 sum_row=0;
453 main_row=0;
454
455 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
456 j=A->pattern->index[iptr];
457 fnorm=0;
458 #pragma ivdep
459 for(bi=0;bi<n_block*n_block;++bi) fnorm+=A->val[iptr*n_block*n_block+bi]*A->val[iptr*n_block*n_block+bi];
460 fnorm=sqrt(fnorm);
461 rtmp[iptr-A->pattern->ptr[i]]=fnorm;
462 if( j != i) {
463 max_offdiagonal = MAX(max_offdiagonal,fnorm);
464 sum_row+=fnorm;
465 } else {
466 main_row=fnorm;
467 }
468 }
469 threshold = theta*max_offdiagonal;
470
471 kdeg=0;
472 if (tau*main_row < sum_row) { /* no diagonal domainance */
473 #pragma ivdep
474 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
475 j=A->pattern->index[iptr];
476 if(rtmp[iptr-A->pattern->ptr[i]] > threshold && i!=j) {
477 S[A->pattern->ptr[i]+kdeg] = j;
478 kdeg++;
479 }
480 }
481 }
482 degree_S[i]=kdeg;
483 }
484 TMPMEMFREE(rtmp);
485 } /* end of parallel region */
486
487 }
488
489 /* the runge stueben coarsening algorithm: */
490 void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset_S,
491 const dim_t* degree_S, const index_t* S,
492 index_t*split_marker, const bool_t usePanel)
493 {
494
495 index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;
496 dim_t i,k, p, q, *degree_ST=NULL, len_panel, len_panel_new;
497 register index_t j, itmp;
498
499 if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
500
501 lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);
502 degree_ST=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degree_ST);
503 ST=TMPMEMALLOC(offset_S[n], index_t); Esys_checkPtr(ST);
504 if (usePanel) {
505 notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);
506 panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);
507 }
508
509
510
511 if (Esys_noError() ) {
512 /* initialize split_marker and split_marker :*/
513 /* those unknows which are not influenced go into F, the rest is available for F or C */
514 #pragma omp parallel for private(i) schedule(static)
515 for (i=0;i<n;++i) {
516 degree_ST[i]=0;
517 if (degree_S[i]>0) {
518 lambda[i]=0;
519 split_marker[i]=PASO_AMG_UNDECIDED;
520 } else {
521 split_marker[i]=PASO_AMG_IN_F;
522 lambda[i]=-1;
523 }
524 }
525 /* create transpose :*/
526 for (i=0;i<n;++i) {
527 for (p=0; p<degree_S[i]; ++p) {
528 j=S[offset_S[i]+p];
529 ST[offset_S[j]+degree_ST[j]]=i;
530 degree_ST[j]++;
531 }
532 }
533 /* lambda[i] = |undecided k in ST[i]| + 2 * |F-unknown in ST[i]| */
534 #pragma omp parallel for private(i, j, itmp) schedule(static)
535 for (i=0;i<n;++i) {
536 if (split_marker[i]==PASO_AMG_UNDECIDED) {
537 itmp=lambda[i];
538 for (p=0; p<degree_ST[i]; ++p) {
539 j=ST[offset_S[i]+p];
540 if (split_marker[j]==PASO_AMG_UNDECIDED) {
541 itmp++;
542 } else { /* at this point there are no C points */
543 itmp+=2;
544 }
545 }
546 lambda[i]=itmp;
547 }
548 }
549 if (usePanel) {
550 #pragma omp parallel for private(i) schedule(static)
551 for (i=0;i<n;++i) notInPanel[i]=TRUE;
552 }
553 /* start search :*/
554 i=Paso_Util_arg_max(n,lambda);
555 while (lambda[i]>-1) { /* is there any undecided unknown? */
556
557 if (usePanel) {
558 len_panel=0;
559 do {
560 /* the unknown i is moved to C */
561 split_marker[i]=PASO_AMG_IN_C;
562 lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
563
564 /* all undecided unknown strongly coupled to i are moved to F */
565 for (p=0; p<degree_ST[i]; ++p) {
566 j=ST[offset_S[i]+p];
567
568 if (split_marker[j]==PASO_AMG_UNDECIDED) {
569
570 split_marker[j]=PASO_AMG_IN_F;
571 lambda[j]=-1;
572
573 for (q=0; q<degree_ST[j]; ++q) {
574 k=ST[offset_S[j]+q];
575 if (split_marker[k]==PASO_AMG_UNDECIDED) {
576 lambda[k]++;
577 if (notInPanel[k]) {
578 notInPanel[k]=FALSE;
579 panel[len_panel]=k;
580 len_panel++;
581 }
582
583 } /* the unknown i is moved to C */
584 split_marker[i]=PASO_AMG_IN_C;
585 lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
586 }
587
588 }
589 }
590 for (p=0; p<degree_S[i]; ++p) {
591 j=S[offset_S[i]+p];
592 if (split_marker[j]==PASO_AMG_UNDECIDED) {
593 lambda[j]--;
594 if (notInPanel[j]) {
595 notInPanel[j]=FALSE;
596 panel[len_panel]=j;
597 len_panel++;
598 }
599 }
600 }
601
602 /* consolidate panel */
603 /* remove lambda[q]=-1 */
604 lambda_max=-1;
605 i=-1;
606 len_panel_new=0;
607 for (q=0; q<len_panel; q++) {
608 k=panel[q];
609 lambda_k=lambda[k];
610 if (split_marker[k]==PASO_AMG_UNDECIDED) {
611 panel[len_panel_new]=k;
612 len_panel_new++;
613
614 if (lambda_max == lambda_k) {
615 if (k<i) i=k;
616 } else if (lambda_max < lambda_k) {
617 lambda_max =lambda_k;
618 i=k;
619 }
620 }
621 }
622 len_panel=len_panel_new;
623 } while (len_panel>0);
624 } else {
625 /* the unknown i is moved to C */
626 split_marker[i]=PASO_AMG_IN_C;
627 lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
628
629 /* all undecided unknown strongly coupled to i are moved to F */
630 for (p=0; p<degree_ST[i]; ++p) {
631 j=ST[offset_S[i]+p];
632 if (split_marker[j]==PASO_AMG_UNDECIDED) {
633
634 split_marker[j]=PASO_AMG_IN_F;
635 lambda[j]=-1;
636
637 for (q=0; q<degree_ST[j]; ++q) {
638 k=ST[offset_S[j]+q];
639 if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;
640 }
641
642 }
643 }
644 for (p=0; p<degree_S[i]; ++p) {
645 j=S[offset_S[i]+p];
646 if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;
647 }
648
649 }
650 i=Paso_Util_arg_max(n,lambda);
651 }
652
653 }
654 TMPMEMFREE(lambda);
655 TMPMEMFREE(ST);
656 TMPMEMFREE(degree_ST);
657 TMPMEMFREE(panel);
658 TMPMEMFREE(notInPanel);
659 }
660 /* ensures that two F nodes are connected via a C node :*/
661 void Paso_Preconditioner_AMG_enforceFFConnectivity(const dim_t n, const index_t* offset_S,
662 const dim_t* degree_S, const index_t* S,
663 index_t*split_marker)
664 {
665 dim_t i, p, q;
666
667 /* now we make sure that two (strongly) connected F nodes are (strongly) connected via a C node. */
668 for (i=0;i<n;++i) {
669 if ( (split_marker[i]==PASO_AMG_IN_F) && (degree_S[i]>0) ) {
670 for (p=0; p<degree_S[i]; ++p) {
671 register index_t j=S[offset_S[i]+p];
672 if ( (split_marker[j]==PASO_AMG_IN_F) && (degree_S[j]>0) ) {
673 /* i and j are now two F nodes which are strongly connected */
674 /* is there a C node they share ? */
675 register index_t sharing=-1;
676 for (q=0; q<degree_S[i]; ++q) {
677 index_t k=S[offset_S[i]+q];
678 if (split_marker[k]==PASO_AMG_IN_C) {
679 register index_t* where_k=(index_t*)bsearch(&k, &(S[offset_S[j]]), degree_S[j], sizeof(index_t), Paso_comparIndex);
680 if (where_k != NULL) {
681 sharing=k;
682 break;
683 }
684 }
685 }
686 if (sharing<0) {
687 if (i<j) {
688 split_marker[j]=PASO_AMG_IN_C;
689 } else {
690 split_marker[i]=PASO_AMG_IN_C;
691 break; /* no point to look any further as i is now a C node */
692 }
693 }
694 }
695 }
696 }
697 }
698 }

  ViewVC Help
Powered by ViewVC 1.1.26