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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26