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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3403 - (show 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
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* split_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F=NULL, *S=NULL, *degree=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 if (verbose) printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
109 level, options->level_max, A_p->pattern->len/(1.*n * n), options->min_coarse_sparsity, n, options->min_coarse_matrix_size );
110 return NULL;
111 }
112 /* Start Coarsening : */
113
114 split_marker=TMPMEMALLOC(n,index_t);
115 counter=TMPMEMALLOC(n,index_t);
116 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 /*
120 set splitting of unknows:
121
122 */
123 time0=Esys_timer();
124 if (n_block>1) {
125 Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree, S, theta,tau);
126 } else {
127 Paso_Preconditioner_AMG_setStrongConnections(A_p, degree, S, theta,tau);
128 }
129 Paso_Preconditioner_AMG_RungeStuebenSearch(n, A_p->pattern->ptr, degree, S, split_marker, options->usePanel);
130 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 if (verbose) printf("Paso_Preconditioner: AMG level %d: %d unknowns are flagged for elimination. %d left.\n",level,n_F,n-n_F);
142
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 if (! Esys_checkPtr(out)) {
148 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 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 out->Smoother = Paso_Preconditioner_LocalSmoother_alloc(A_p, (options->smoother == PASO_JACOBI), verbose);
170
171 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
174 /* 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 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
195 #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 }
203 /*
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 }
210 /*
211 construct Prolongation operator as transposed of restriction operator:
212 */
213 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 /*
219 construct coarse level matrix:
220 */
221 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
230 /*
231 constructe courser level:
232
233 */
234 if ( Esys_noError()) {
235 out->AMG_C=Paso_Preconditioner_LocalAMG_alloc(A_C,level+1,options);
236 }
237 if ( Esys_noError()) {
238 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 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_C);
244 Paso_SparseMatrix_free(A_C);
245 out->A_C->solver_package = PASO_MKL;
246 if (verbose) printf("Paso_Preconditioner: AMG: use MKL direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
247 #else
248 #ifdef UMFPACK
249 out->A_C=Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_C);
250 Paso_SparseMatrix_free(A_C);
251 out->A_C->solver_package = PASO_UMFPACK;
252 if (verbose) printf("Paso_Preconditioner: AMG: use UMFPACK direct solver on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
253 #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 if (verbose) printf("Paso_Preconditioner: AMG: use smoother on the coarsest level (number of unknowns = %d).\n",n_C*n_block);
258 #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 }
272 TMPMEMFREE(counter);
273 TMPMEMFREE(split_marker);
274 TMPMEMFREE(degree);
275 TMPMEMFREE(S);
276
277 if (Esys_noError()) {
278 return out;
279 } else {
280 Paso_Preconditioner_LocalAMG_free(out);
281 return NULL;
282 }
283 }
284
285
286 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
292 /* 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 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 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 Paso_Preconditioner_LocalSmoother_solve(amg->A_C, amg->A_C->solver_p,amg->x_C,amg->b_C,pre_sweeps+post_sweeps, FALSE);
318 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 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */
326
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
340 /* theta = threshold for strong connections */
341 /* tau = threshold for diagonal dominance */
342
343 /*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 {
352 const dim_t n=A->numRows;
353 index_t iptr, i,j;
354 dim_t kdeg;
355 double max_offdiagonal, threshold, sum_row, main_row, fnorm;
356
357
358 #pragma omp parallel for private(i,iptr,max_offdiagonal, threshold,j, kdeg, sum_row, main_row, fnorm) schedule(static)
359 for (i=0;i<n;++i) {
360 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
378 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
391 }
392
393 /* theta = threshold for strong connections */
394 /* tau = threshold for diagonal dominance */
395 /*S_i={j \in N_i; i strongly coupled to j}
396
397 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 {
404 const dim_t n_block=A->row_block_size;
405 const dim_t n=A->numRows;
406 index_t iptr, i,j, bi;
407 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
418 rtmp=TMPMEMALLOC(max_deg, double);
419
420 #pragma omp for schedule(static)
421 for (i=0;i<n;++i) {
422
423 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 rtmp[iptr-A->pattern->ptr[i]]=fnorm;
433 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
442 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
458 }
459
460 /* the runge stueben coarsening algorithm: */
461 void Paso_Preconditioner_AMG_RungeStuebenSearch(const dim_t n, const index_t* offset,
462 const dim_t* degree, const index_t* S,
463 index_t*split_marker, const bool_t usePanel)
464 {
465
466 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 if (n<=0) return; /* make sure that the return of Paso_Util_arg_max is not pointing to nirvana */
471
472 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
480
481
482 if (Esys_noError() ) {
483 /* 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 if (usePanel) {
521 #pragma omp parallel for private(i) schedule(static)
522 for (i=0;i<n;++i) notInPanel[i]=TRUE;
523 }
524 /* start search :*/
525 i=Paso_Util_arg_max(n,lambda);
526 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 }
561 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
573 /* 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 }
615 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 }
621 i=Paso_Util_arg_max(n,lambda);
622 }
623 }
624 TMPMEMFREE(lambda);
625 TMPMEMFREE(ST);
626 TMPMEMFREE(degreeT);
627 TMPMEMFREE(panel);
628 TMPMEMFREE(notInPanel);
629 }

  ViewVC Help
Powered by ViewVC 1.1.26