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

Annotation of /trunk/paso/src/LocalAMG.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3441 - (hide annotations)
Fri Jan 14 01:09:09 2011 UTC (9 years, 10 months ago) by gross
File MIME type: text/plain
File size: 22806 byte(s)
clarifiaction of function names as a first step towards a MPI version of AMG
1 gross 3441
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 (local version) */
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_LocalAMG_setStrongConnections_Block(A_p, degree_S, S, theta,tau);
150     } else {
151     Paso_Preconditioner_LocalAMG_setStrongConnections(A_p, degree_S, S, theta,tau);
152     }
153     Paso_Preconditioner_LocalAMG_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_LocalAMG_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_LocalAMG_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_LocalAMG_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_LocalAMG_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_LocalAMG_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_LocalAMG_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