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

Contents of /trunk/paso/src/LocalAMG.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4803 - (show annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 10 months ago) by caltinay
File size: 24979 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

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

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/LocalAMG.cpp:3531-3826 /branches/lapack2681/paso/src/LocalAMG.cpp:2682-2741 /branches/pasowrap/paso/src/LocalAMG.cpp:3661-3674 /branches/py3_attempt2/paso/src/LocalAMG.cpp:3871-3891 /branches/restext/paso/src/LocalAMG.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/LocalAMG.cpp:3669-3791 /branches/stage3.0/paso/src/LocalAMG.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/LocalAMG.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/LocalAMG.cpp:3517-3974 /release/3.0/paso/src/LocalAMG.cpp:2591-2601 /trunk/paso/src/LocalAMG.cpp:4257-4344 /trunk/ripley/test/python/paso/src/LocalAMG.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26