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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3402 - (hide annotations)
Tue Dec 7 07:36:12 2010 UTC (8 years, 10 months ago) by gross
File MIME type: text/plain
File size: 19830 byte(s)
AMG support now block size > 3 (requires clapack or MKL).

additional diagnostics for AMG 


1 gross 3193
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 gross 3303 /* Author: artak@uq.edu.au, l.gross@uq.edu.au */
22 gross 3193
23     /**************************************************************/
24    
25 gross 3283 #define SHOW_TIMING FALSE
26    
27 gross 3193 #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 gross 3283 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 gross 3193
49 gross 3283
50     MEMFREE(in);
51 gross 3193 }
52     }
53    
54 gross 3402 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     return DBLE(in->A_C->pattern->len)/DBLE(in->A_C->numRows)/DBLE(in->A_C->numRows);
64     } else {
65     return Paso_Preconditioner_LocalAMG_getCoarseLevelSparsity(in->AMG_C);
66     }
67     }
68     dim_t Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(const Paso_Preconditioner_LocalAMG * in) {
69     if (in->AMG_C == NULL) {
70     return in->A_C->numRows;
71     } else {
72     return Paso_Preconditioner_LocalAMG_getNumCoarseUnknwons(in->AMG_C);
73     }
74     }
75 gross 3283 /*****************************************************************
76 gross 3193
77 gross 3283 constructs AMG
78    
79     ******************************************************************/
80     Paso_Preconditioner_LocalAMG* Paso_Preconditioner_LocalAMG_alloc(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
81 gross 3193
82     Paso_Preconditioner_LocalAMG* out=NULL;
83     bool_t verbose=options->verbose;
84    
85 gross 3303 Paso_SparseMatrix *Atemp=NULL, *A_C=NULL;
86 gross 3193 const dim_t n=A_p->numRows;
87     const dim_t n_block=A_p->row_block_size;
88 gross 3303 index_t* split_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree=NULL;
89 gross 3283 dim_t n_F=0, n_C=0, i;
90     double time0=0;
91     const double theta = options->coarsening_threshold;
92     const double tau = options->diagonal_dominance_threshold;
93 gross 3193
94    
95     /*
96 gross 3283 is the input matrix A suitable for coarsening
97    
98 gross 3193 */
99 gross 3283 if ( (A_p->pattern->len >= options->min_coarse_sparsity * n * n ) || (n <= options->min_coarse_matrix_size) || (level > options->level_max) ) {
100 gross 3402 if (verbose) printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
101 gross 3283 level, options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size );
102 gross 3193 return NULL;
103 gross 3283 }
104     /* Start Coarsening : */
105 gross 3193
106 gross 3283 split_marker=TMPMEMALLOC(n,index_t);
107     counter=TMPMEMALLOC(n,index_t);
108 gross 3303 degree=TMPMEMALLOC(n, dim_t);
109     S=TMPMEMALLOC(A_p->pattern->len, index_t);
110     if ( !( Esys_checkPtr(split_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree) || Esys_checkPtr(S) ) ) {
111 gross 3283 /*
112     set splitting of unknows:
113    
114     */
115     time0=Esys_timer();
116     if (n_block>1) {
117 gross 3303 Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree, S, theta,tau);
118 gross 3283 } else {
119 gross 3303 Paso_Preconditioner_AMG_setStrongConnections(A_p, degree, S, theta,tau);
120 gross 3283 }
121 gross 3352 Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree, S, split_marker, options->usePanel);
122 gross 3283 options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
123    
124     if (Esys_noError() ) {
125     #pragma omp parallel for private(i) schedule(static)
126     for (i = 0; i < n; ++i) split_marker[i]= (split_marker[i] == PASO_AMG_IN_F);
127    
128     /*
129     count number of unkowns to be eliminated:
130     */
131     n_F=Paso_Util_cumsum_maskedTrue(n,counter, split_marker);
132     n_C=n-n_F;
133 gross 3402 if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
134 gross 3283
135     if ( n_F == 0 ) { /* is a nasty case. a direct solver should be used, return NULL */
136     out = NULL;
137     } else {
138     out=MEMALLOC(1,Paso_Preconditioner_LocalAMG);
139 gross 3323 if (! Esys_checkPtr(out)) {
140 gross 3283 out->level = level;
141     out->n = n;
142     out->n_F = n_F;
143     out->n_block = n_block;
144     out->A_C = NULL;
145     out->P = NULL;
146     out->R = NULL;
147     out->post_sweeps = options->post_sweeps;
148     out->pre_sweeps = options->pre_sweeps;
149     out->r = NULL;
150     out->x_C = NULL;
151     out->b_C = NULL;
152     out->AMG_C = NULL;
153 gross 3323 out->Smoother=NULL;
154     }
155     mask_C=TMPMEMALLOC(n,index_t);
156     rows_in_F=TMPMEMALLOC(n_F,index_t);
157     Esys_checkPtr(mask_C);
158     Esys_checkPtr(rows_in_F);
159     if ( Esys_noError() ) {
160    
161 gross 3283 out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
162 gross 3323
163 gross 3402 if (n_C != 0) {
164     /* if nothing is been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
165 gross 3283
166 gross 3315 /* allocate helpers :*/
167     out->x_C=MEMALLOC(n_block*n_C,double);
168     out->b_C=MEMALLOC(n_block*n_C,double);
169     out->r=MEMALLOC(n_block*n,double);
170    
171     Esys_checkPtr(out->r);
172     Esys_checkPtr(out->x_C);
173     Esys_checkPtr(out->b_C);
174    
175 gross 3323 if ( Esys_noError() ) {
176     /* creates index for F:*/
177     #pragma omp parallel private(i)
178     {
179     #pragma omp for schedule(static)
180     for (i = 0; i < n; ++i) {
181     if (split_marker[i]) rows_in_F[counter[i]]=i;
182     }
183     }
184     /* create mask of C nodes with value >-1 gives new id */
185     i=Paso_Util_cumsum_maskedFalse(n,counter, split_marker);
186 gross 3193
187 gross 3323 #pragma omp parallel for private(i) schedule(static)
188     for (i = 0; i < n; ++i) {
189     if (split_marker[i]) {
190     mask_C[i]=-1;
191     } else {
192     mask_C[i]=counter[i];;
193     }
194 gross 3283 }
195 gross 3323 /*
196     get Restriction :
197     */
198     time0=Esys_timer();
199     out->P=Paso_Preconditioner_AMG_getDirectProlongation(A_p,degree,S,n_C,mask_C);
200     if (SHOW_TIMING) printf("timing: level %d: getProlongation: %e\n",level, Esys_timer()-time0);
201 gross 3283 }
202 gross 3323 /*
203 gross 3283 construct Prolongation operator as transposed of restriction operator:
204     */
205 gross 3315 if ( Esys_noError()) {
206     time0=Esys_timer();
207     out->R=Paso_SparseMatrix_getTranspose(out->P);
208     if (SHOW_TIMING) printf("timing: level %d: Paso_SparseMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
209     }
210 gross 3283 /*
211 gross 3312 construct coarse level matrix:
212 gross 3283 */
213 gross 3315 if ( Esys_noError()) {
214     time0=Esys_timer();
215     Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
216     A_C=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
217     Paso_SparseMatrix_free(Atemp);
218     if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);
219     }
220    
221 gross 3283
222     /*
223     constructe courser level:
224    
225     */
226     if ( Esys_noError()) {
227 gross 3315 out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);
228     }
229     if ( Esys_noError()) {
230 gross 3283 if ( out->AMG_C == NULL ) {
231     out->reordering = options->reordering;
232     out->refinements = options->coarse_matrix_refinements;
233     /* no coarse level matrix has been constructed. use direct solver */
234     #ifdef MKL
235 gross 3312 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
236 gross 3283 Paso_SparseMatrix_free(A_C);
237     out->A_C->solver_package = PASO_MKL;
238 gross 3402 if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
239 gross 3283 #else
240     #ifdef UMFPACK
241 gross 3312 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
242 gross 3283 Paso_SparseMatrix_free(A_C);
243     out->A_C->solver_package = PASO_UMFPACK;
244 gross 3402 if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
245 gross 3283 #else
246     out->A_C=A_C;
247     out->A_C->solver_p=Paso_Preconditioner_LocalSmoother_alloc(out->A_C, (options->smoother == PASO_JACOBI), verbose);
248     out->A_C->solver_package = PASO_SMOOTHER;
249 gross 3402 if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
250 gross 3283 #endif
251     #endif
252     } else {
253     /* finally we set some helpers for the solver step */
254     out->A_C=A_C;
255     }
256     }
257     }
258     }
259     TMPMEMFREE(mask_C);
260     TMPMEMFREE(rows_in_F);
261     }
262     }
263 gross 3193 }
264     TMPMEMFREE(counter);
265 gross 3283 TMPMEMFREE(split_marker);
266 gross 3303 TMPMEMFREE(degree);
267     TMPMEMFREE(S);
268 gross 3193
269 jfenwick 3259 if (Esys_noError()) {
270 gross 3193 return out;
271     } else {
272     Paso_Preconditioner_LocalAMG_free(out);
273     return NULL;
274     }
275     }
276    
277    
278 gross 3283 void Paso_Preconditioner_LocalAMG_solve(Paso_SparseMatrix* A, Paso_Preconditioner_LocalAMG * amg, double * x, double * b) {
279     const dim_t n = amg->n * amg->n_block;
280     double time0=0;
281     const dim_t post_sweeps=amg->post_sweeps;
282     const dim_t pre_sweeps=amg->pre_sweeps;
283 gross 3193
284 gross 3283 /* presmoothing */
285     time0=Esys_timer();
286     Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
287     time0=Esys_timer()-time0;
288     if (SHOW_TIMING) printf("timing: level %d: Presmooting: %e\n",amg->level, time0);
289     /* end of presmoothing */
290    
291     if (amg->n_F < amg->n) { /* is there work on the coarse level? */
292     time0=Esys_timer();
293     Paso_Copy(n, amg->r, b); /* r <- b */
294 gross 3312 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
295     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->R,amg->r,0.,amg->b_C); /* b_c = R*r */
296 gross 3283 time0=Esys_timer()-time0;
297     /* coarse level solve */
298     if ( amg->AMG_C == NULL) {
299     time0=Esys_timer();
300     /* A_C is the coarsest level */
301     switch (amg->A_C->solver_package) {
302     case (PASO_MKL):
303     Paso_MKL(amg->A_C, amg->x_C,amg->b_C, amg->reordering, amg->refinements, SHOW_TIMING);
304     break;
305     case (PASO_UMFPACK):
306     Paso_UMFPACK(amg->A_C, amg->x_C,amg->b_C, amg->refinements, SHOW_TIMING);
307     break;
308     case (PASO_SMOOTHER):
309 gross 3402 Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);
310 gross 3283 break;
311     }
312     if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
313     } else {
314     Paso_Preconditioner_LocalAMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C) */
315     }
316     time0=time0+Esys_timer();
317 gross 3303 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */
318 gross 3283
319     /*postsmoothing*/
320    
321     /*solve Ax=b with initial guess x */
322     time0=Esys_timer();
323     Paso_Preconditioner_LocalSmoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
324     time0=Esys_timer()-time0;
325     if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
326     /*end of postsmoothing*/
327    
328     }
329     return;
330     }
331 gross 3193
332 gross 3283 /* theta = threshold for strong connections */
333     /* tau = threshold for diagonal dominance */
334 gross 3193
335 gross 3303 /*S_i={j \in N_i; i strongly coupled to j}
336    
337     in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
338     */
339    
340     void Paso_Preconditioner_AMG_setStrongConnections(Paso_SparseMatrix* A,
341     dim_t *degree, index_t *S,
342     const double theta, const double tau)
343 gross 3283 {
344     const dim_t n=A->numRows;
345 gross 3303 index_t iptr, i,j;
346 gross 3283 dim_t kdeg;
347     double max_offdiagonal, threshold, sum_row, main_row, fnorm;
348 gross 3193
349 gross 3303
350 gross 3283 #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)
351 gross 3402 for (i=0;i<n;++i) {
352 gross 3283 max_offdiagonal = 0.;
353     sum_row=0;
354     main_row=0;
355     #pragma ivdep
356     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
357     j=A->pattern->index[iptr];
358     fnorm=ABS(A->val[iptr]);
359    
360     if( j != i) {
361     max_offdiagonal = MAX(max_offdiagonal,fnorm);
362     sum_row+=fnorm;
363     } else {
364     main_row=fnorm;
365     }
366     }
367     threshold = theta*max_offdiagonal;
368     kdeg=0;
369 gross 3402
370 gross 3283 if (tau*main_row < sum_row) { /* no diagonal domainance */
371     #pragma ivdep
372     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
373     j=A->pattern->index[iptr];
374     if(ABS(A->val[iptr])>threshold && i!=j) {
375     S[A->pattern->ptr[i]+kdeg] = j;
376     kdeg++;
377     }
378     }
379     }
380     degree[i]=kdeg;
381     }
382 gross 3303
383 gross 3283 }
384 gross 3193
385 gross 3283 /* theta = threshold for strong connections */
386     /* tau = threshold for diagonal dominance */
387 gross 3303 /*S_i={j \in N_i; i strongly coupled to j}
388 gross 3283
389 gross 3303 in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
390     */
391     void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SparseMatrix* A,
392     dim_t *degree, index_t *S,
393     const double theta, const double tau)
394    
395 gross 3283 {
396     const dim_t n_block=A->row_block_size;
397     const dim_t n=A->numRows;
398 gross 3303 index_t iptr, i,j, bi;
399 gross 3283 dim_t kdeg, max_deg;
400     register double max_offdiagonal, threshold, fnorm, sum_row, main_row;
401     double *rtmp;
402    
403    
404     #pragma omp parallel private(i,iptr,max_offdiagonal, kdeg, threshold,j, max_deg, fnorm, sum_row, main_row, rtmp)
405     {
406     max_deg=0;
407     #pragma omp for schedule(static)
408     for (i=0;i<n;++i) max_deg=MAX(max_deg, A->pattern->ptr[i+1]-A->pattern->ptr[i]);
409 gross 3193
410 gross 3283 rtmp=TMPMEMALLOC(max_deg, double);
411    
412     #pragma omp for schedule(static)
413     for (i=0;i<n;++i) {
414 gross 3193
415 gross 3283 max_offdiagonal = 0.;
416     sum_row=0;
417     main_row=0;
418     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
419     j=A->pattern->index[iptr];
420     fnorm=0;
421     #pragma ivdep
422     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];
423     fnorm=sqrt(fnorm);
424 gross 3287 rtmp[iptr-A->pattern->ptr[i]]=fnorm;
425 gross 3283 if( j != i) {
426     max_offdiagonal = MAX(max_offdiagonal,fnorm);
427     sum_row+=fnorm;
428     } else {
429     main_row=fnorm;
430     }
431     }
432     threshold = theta*max_offdiagonal;
433 gross 3193
434 gross 3283 kdeg=0;
435     if (tau*main_row < sum_row) { /* no diagonal domainance */
436     #pragma ivdep
437     for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
438     j=A->pattern->index[iptr];
439     if(rtmp[iptr-A->pattern->ptr[i]] > threshold && i!=j) {
440     S[A->pattern->ptr[i]+kdeg] = j;
441     kdeg++;
442     }
443     }
444     }
445     degree[i]=kdeg;
446     }
447     TMPMEMFREE(rtmp);
448     } /* end of parallel region */
449 gross 3303
450 gross 3283 }
451    
452     /* the runge stueben coarsening algorithm: */
453 gross 3303 void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset,
454 gross 3352 const dim_t* degree, const index_t* S,
455     index_t*split_marker, const bool_t usePanel)
456 gross 3283 {
457    
458 gross 3351 index_t *lambda=NULL, *ST=NULL, *notInPanel=NULL, *panel=NULL, lambda_max, lambda_k;
459     dim_t i,k, p, q, *degreeT=NULL, len_panel, len_panel_new;
460     register index_t j, itmp;
461    
462 gross 3283 if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
463    
464 gross 3351 lambda=TMPMEMALLOC(n, index_t); Esys_checkPtr(lambda);
465     degreeT=TMPMEMALLOC(n, dim_t); Esys_checkPtr(degreeT);
466     ST=TMPMEMALLOC(offset[n], index_t); Esys_checkPtr(ST);
467     if (usePanel) {
468     notInPanel=TMPMEMALLOC(n, bool_t); Esys_checkPtr(notInPanel);
469     panel=TMPMEMALLOC(n, index_t); Esys_checkPtr(panel);
470     }
471 gross 3283
472 gross 3351
473    
474     if (Esys_noError() ) {
475 gross 3283 /* initialize split_marker and split_marker :*/
476     /* those unknows which are not influenced go into F, the rest is available for F or C */
477     #pragma omp parallel for private(i) schedule(static)
478     for (i=0;i<n;++i) {
479     degreeT[i]=0;
480     if (degree[i]>0) {
481     lambda[i]=0;
482     split_marker[i]=PASO_AMG_UNDECIDED;
483     } else {
484     split_marker[i]=PASO_AMG_IN_F;
485     lambda[i]=-1;
486     }
487     }
488     /* create transpose :*/
489     for (i=0;i<n;++i) {
490     for (p=0; p<degree[i]; ++p) {
491     j=S[offset[i]+p];
492     ST[offset[j]+degreeT[j]]=i;
493     degreeT[j]++;
494     }
495     }
496     /* lambda[i] = |undecided k in ST[i]| + 2 * |F-unknown in ST[i]| */
497     #pragma omp parallel for private(i, j, itmp) schedule(static)
498     for (i=0;i<n;++i) {
499     if (split_marker[i]==PASO_AMG_UNDECIDED) {
500     itmp=lambda[i];
501     for (p=0; p<degreeT[i]; ++p) {
502     j=ST[offset[i]+p];
503     if (split_marker[j]==PASO_AMG_UNDECIDED) {
504     itmp++;
505     } else { /* at this point there are no C points */
506     itmp+=2;
507     }
508     }
509     lambda[i]=itmp;
510     }
511     }
512 gross 3351 if (usePanel) {
513     #pragma omp parallel for private(i) schedule(static)
514     for (i=0;i<n;++i) notInPanel[i]=TRUE;
515     }
516 gross 3283 /* start search :*/
517     i=Paso_Util_arg_max(n,lambda);
518 gross 3351 while (lambda[i]>-1) { /* is there any undecided unknown? */
519    
520     if (usePanel) {
521     len_panel=0;
522     do {
523     /* the unknown i is moved to C */
524     split_marker[i]=PASO_AMG_IN_C;
525     lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
526    
527     /* all undecided unknown strongly coupled to i are moved to F */
528     for (p=0; p<degreeT[i]; ++p) {
529     j=ST[offset[i]+p];
530    
531     if (split_marker[j]==PASO_AMG_UNDECIDED) {
532    
533     split_marker[j]=PASO_AMG_IN_F;
534     lambda[j]=-1;
535    
536     for (q=0; q<degreeT[j]; ++q) {
537     k=ST[offset[j]+q];
538     if (split_marker[k]==PASO_AMG_UNDECIDED) {
539     lambda[k]++;
540     if (notInPanel[k]) {
541     notInPanel[k]=FALSE;
542     panel[len_panel]=k;
543     len_panel++;
544     }
545    
546     } /* the unknown i is moved to C */
547     split_marker[i]=PASO_AMG_IN_C;
548     lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
549     }
550    
551     }
552 gross 3283 }
553 gross 3351 for (p=0; p<degree[i]; ++p) {
554     j=S[offset[i]+p];
555     if (split_marker[j]==PASO_AMG_UNDECIDED) {
556     lambda[j]--;
557     if (notInPanel[j]) {
558     notInPanel[j]=FALSE;
559     panel[len_panel]=j;
560     len_panel++;
561     }
562     }
563     }
564 gross 3193
565 gross 3351 /* consolidate panel */
566     /* remove lambda[q]=-1 */
567     lambda_max=-1;
568     i=-1;
569     len_panel_new=0;
570     for (q=0; q<len_panel; q++) {
571     k=panel[q];
572     lambda_k=lambda[k];
573     if (split_marker[k]==PASO_AMG_UNDECIDED) {
574     panel[len_panel_new]=k;
575     len_panel_new++;
576    
577     if (lambda_max == lambda_k) {
578     if (k<i) i=k;
579     } else if (lambda_max < lambda_k) {
580     lambda_max =lambda_k;
581     i=k;
582     }
583     }
584     }
585     len_panel=len_panel_new;
586     } while (len_panel>0);
587     } else {
588     /* the unknown i is moved to C */
589     split_marker[i]=PASO_AMG_IN_C;
590     lambda[i]=-1; /* lambda from unavailable unknowns is set to -1 */
591    
592     /* all undecided unknown strongly coupled to i are moved to F */
593     for (p=0; p<degreeT[i]; ++p) {
594     j=ST[offset[i]+p];
595     if (split_marker[j]==PASO_AMG_UNDECIDED) {
596    
597     split_marker[j]=PASO_AMG_IN_F;
598     lambda[j]=-1;
599    
600     for (q=0; q<degreeT[j]; ++q) {
601     k=ST[offset[j]+q];
602     if (split_marker[k]==PASO_AMG_UNDECIDED) lambda[k]++;
603     }
604    
605     }
606 gross 3283 }
607 gross 3351 for (p=0; p<degree[i]; ++p) {
608     j=S[offset[i]+p];
609     if(split_marker[j]==PASO_AMG_UNDECIDED) lambda[j]--;
610     }
611    
612 gross 3283 }
613     i=Paso_Util_arg_max(n,lambda);
614     }
615     }
616     TMPMEMFREE(lambda);
617     TMPMEMFREE(ST);
618     TMPMEMFREE(degreeT);
619 gross 3351 TMPMEMFREE(panel);
620     TMPMEMFREE(notInPanel);
621 gross 3193 }

  ViewVC Help
Powered by ViewVC 1.1.26