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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3885 - (show annotations)
Wed Apr 4 22:12:26 2012 UTC (7 years, 5 months ago) by gross
File MIME type: text/plain
File size: 37502 byte(s)
some fix in AMG
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 (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 #define MY_DEBUG 0
27 #define MY_DEBUG1 1
28
29 #include "Paso.h"
30 #include "Preconditioner.h"
31 #include "Options.h"
32 #include "PasoUtil.h"
33 #include "UMFPACK.h"
34 #include "MKL.h"
35 #include<stdio.h>
36
37
38 /**************************************************************/
39
40 /* free all memory used by AMG */
41
42 void Paso_Preconditioner_AMG_free(Paso_Preconditioner_AMG * in) {
43 if (in!=NULL) {
44 Paso_Preconditioner_Smoother_free(in->Smoother);
45 Paso_SystemMatrix_free(in->P);
46 Paso_SystemMatrix_free(in->R);
47 Paso_SystemMatrix_free(in->A_C);
48 Paso_Preconditioner_AMG_free(in->AMG_C);
49 MEMFREE(in->r);
50 MEMFREE(in->x_C);
51 MEMFREE(in->b_C);
52
53 MEMFREE(in);
54 }
55 }
56
57 index_t Paso_Preconditioner_AMG_getMaxLevel(const Paso_Preconditioner_AMG * in) {
58 if (in->AMG_C == NULL) {
59 return in->level;
60 } else {
61 return Paso_Preconditioner_AMG_getMaxLevel(in->AMG_C);
62 }
63 }
64 double Paso_Preconditioner_AMG_getCoarseLevelSparsity(const Paso_Preconditioner_AMG * in) {
65 if (in->AMG_C == NULL) {
66 if (in->A_C == NULL) {
67 return 1.;
68 } else {
69 return Paso_SystemMatrix_getSparsity(in->A_C);
70 }
71 } else {
72 return Paso_Preconditioner_AMG_getCoarseLevelSparsity(in->AMG_C);
73 }
74 }
75 dim_t Paso_Preconditioner_AMG_getNumCoarseUnknwons(const Paso_Preconditioner_AMG * in) {
76 if (in->AMG_C == NULL) {
77 if (in->A_C == NULL) {
78 return 0;
79 } else {
80 return Paso_SystemMatrix_getTotalNumRows(in->A_C);
81 }
82 } else {
83 return Paso_Preconditioner_AMG_getNumCoarseUnknwons(in->AMG_C);
84 }
85 }
86 /*****************************************************************
87
88 constructs AMG
89
90 ******************************************************************/
91 Paso_Preconditioner_AMG* Paso_Preconditioner_AMG_alloc(Paso_SystemMatrix *A_p,dim_t level,Paso_Options* options) {
92
93 Paso_Preconditioner_AMG* out=NULL;
94 Paso_SystemMatrix *A_C=NULL;
95 bool_t verbose=options->verbose;
96
97 const dim_t my_n=A_p->mainBlock->numRows;
98 const dim_t overlap_n=A_p->row_coupleBlock->numRows;
99
100 const dim_t n = my_n + overlap_n;
101
102 const dim_t n_block=A_p->row_block_size;
103 index_t* F_marker=NULL, *counter=NULL, *mask_C=NULL, *rows_in_F;
104 dim_t i, n_F, n_C, F_flag, *F_set=NULL, global_n_C=0, global_n_F=0;
105 double time0=0;
106 const double theta = options->coarsening_threshold;
107 const double tau = options->diagonal_dominance_threshold;
108 const double sparsity=Paso_SystemMatrix_getSparsity(A_p);
109 const dim_t global_n=Paso_SystemMatrix_getGlobalNumRows(A_p);
110
111
112 /*
113 is the input matrix A suitable for coarsening?
114
115 */
116 if ( (sparsity >= options->min_coarse_sparsity) ||
117 (global_n <= options->min_coarse_matrix_size) ||
118 (level > options->level_max) ) {
119
120 if (verbose) {
121 /*
122 print stopping condition:
123 - 'SPAR' = min_coarse_matrix_sparsity exceeded
124 - 'SIZE' = min_coarse_matrix_size exceeded
125 - 'LEVEL' = level_max exceeded
126 */
127 printf("Paso_Preconditioner: AMG: termination of coarsening by ");
128
129 if (sparsity >= options->min_coarse_sparsity)
130 printf("SPAR");
131
132 if (global_n <= options->min_coarse_matrix_size)
133 printf("SIZE");
134
135 if (level > options->level_max)
136 printf("LEVEL");
137
138 printf("\n");
139
140 printf("Paso_Preconditioner: AMG level %d (limit = %d) stopped. sparsity = %e (limit = %e), unknowns = %d (limit = %d)\n",
141 level, options->level_max, sparsity, options->min_coarse_sparsity, global_n, options->min_coarse_matrix_size);
142
143 }
144
145 return NULL;
146 } else {
147 /* Start Coarsening : */
148
149 /* this is the table for strong connections combining mainBlock, col_coupleBlock and row_coupleBlock */
150 const dim_t len_S=A_p->mainBlock->pattern->len + A_p->col_coupleBlock->pattern->len + A_p->row_coupleBlock->pattern->len + A_p->row_coupleBlock->numRows * A_p->col_coupleBlock->numCols;
151
152 dim_t* degree_S=TMPMEMALLOC(n, dim_t);
153 index_t *offset_S=TMPMEMALLOC(n, index_t);
154 index_t *S=TMPMEMALLOC(len_S, index_t);
155 dim_t* degree_ST=TMPMEMALLOC(n, dim_t);
156 index_t *offset_ST=TMPMEMALLOC(n, index_t);
157 index_t *ST=TMPMEMALLOC(len_S, index_t);
158
159
160 F_marker=TMPMEMALLOC(n,index_t);
161 counter=TMPMEMALLOC(n,index_t);
162
163 if ( !( Esys_checkPtr(F_marker) || Esys_checkPtr(counter) || Esys_checkPtr(degree_S) || Esys_checkPtr(offset_S) || Esys_checkPtr(S)
164 || Esys_checkPtr(degree_ST) || Esys_checkPtr(offset_ST) || Esys_checkPtr(ST) ) ) {
165 /*
166 make sure that corresponding values in the row_coupleBlock and col_coupleBlock are identical
167 */
168 Paso_SystemMatrix_copyColCoupleBlock(A_p);
169 Paso_SystemMatrix_copyRemoteCoupleBlock(A_p, FALSE);
170
171 /*
172 set splitting of unknows:
173
174 */
175 time0=Esys_timer();
176 if (n_block>1) {
177 Paso_Preconditioner_AMG_setStrongConnections_Block(A_p, degree_S, offset_S, S, theta,tau);
178 } else {
179 Paso_Preconditioner_AMG_setStrongConnections(A_p, degree_S, offset_S, S, theta,tau);
180 }
181 Paso_Preconditioner_AMG_transposeStrongConnections(n, degree_S, offset_S, S, n, degree_ST, offset_ST, ST);
182 /* Paso_SystemMatrix_extendedRowsForST(A_p, degree_ST, offset_ST, ST);
183 */
184
185 Paso_Preconditioner_AMG_CIJPCoarsening(n,my_n,F_marker,
186 degree_S, offset_S, S, degree_ST, offset_ST, ST,
187 A_p->col_coupler->connector,A_p->col_distribution);
188
189
190 /* in BoomerAMG if interpolation is used FF connectivity is required */
191 /*MPI:
192 if (options->interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
193 Paso_Preconditioner_AMG_enforceFFConnectivity(n, A_p->pattern->ptr, degree_S, S, F_marker);
194 */
195
196 options->coarsening_selection_time=Esys_timer()-time0 + MAX(0, options->coarsening_selection_time);
197 if (Esys_noError() ) {
198 #pragma omp parallel for private(i) schedule(static)
199 for (i = 0; i < n; ++i) F_marker[i]=(F_marker[i] == PASO_AMG_IN_F);
200
201 /*
202 count number of unkowns to be eliminated:
203 */
204 n_F=Paso_Util_cumsum_maskedTrue(n,counter, F_marker);
205 /* collect n_F values on all processes, a direct solver should
206 be used if any n_F value is 0 */
207 F_set = TMPMEMALLOC(A_p->mpi_info->size, dim_t);
208 #ifdef ESYS_MPI
209 MPI_Allgather(&n_F, 1, MPI_INT, F_set, 1, MPI_INT, A_p->mpi_info->comm);
210 #endif
211 global_n_F=0;
212 F_flag = 1;
213 for (i=0; i<A_p->mpi_info->size; i++) {
214 global_n_F+=F_set[i];
215 if (F_set[i] == 0) F_flag = 0;
216 printf(" p [%d]=%d\n",i, F_set[i]);
217 }
218 printf(" global n = %d\n",global_n);
219 TMPMEMFREE(F_set);
220
221 n_C=n-n_F;
222 global_n_C=global_n-global_n_F;
223 if (verbose) printf("Paso_Preconditioner: AMG (non-local) level %d: %d unknowns are flagged for elimination. %d left.\n",level,global_n_F,global_n_C);
224
225
226 /* if ( n_F == 0 ) { is a nasty case. a direct solver should be used, return NULL */
227 if (F_flag == 0) {
228 out = NULL;
229 } else {
230 out=MEMALLOC(1,Paso_Preconditioner_AMG);
231 if (! Esys_checkPtr(out)) {
232 out->level = level;
233 out->n = n;
234 out->n_F = n_F;
235 out->n_block = n_block;
236 out->A_C = NULL;
237 out->P = NULL;
238 out->R = NULL;
239 out->post_sweeps = options->post_sweeps;
240 out->pre_sweeps = options->pre_sweeps;
241 out->r = NULL;
242 out->x_C = NULL;
243 out->b_C = NULL;
244 out->AMG_C = NULL;
245 out->Smoother=NULL;
246 }
247 mask_C=TMPMEMALLOC(n,index_t);
248 rows_in_F=TMPMEMALLOC(n_F,index_t);
249 Esys_checkPtr(mask_C);
250 Esys_checkPtr(rows_in_F);
251 if ( Esys_noError() ) {
252
253 out->Smoother = Paso_Preconditioner_Smoother_alloc(A_p, (options->smoother == PASO_JACOBI), 0, verbose);
254
255 if (global_n_C != 0) {
256 /* if nothing has been removed we have a diagonal dominant matrix and we just run a few steps of the smoother */
257
258 /* allocate helpers :*/
259 out->x_C=MEMALLOC(n_block*n_C,double);
260 out->b_C=MEMALLOC(n_block*n_C,double);
261 out->r=MEMALLOC(n_block*n,double);
262
263 Esys_checkPtr(out->r);
264 Esys_checkPtr(out->x_C);
265 Esys_checkPtr(out->b_C);
266
267 if ( Esys_noError() ) {
268 /* creates index for F:*/
269 #pragma omp parallel private(i)
270 {
271 #pragma omp for schedule(static)
272 for (i = 0; i < n; ++i) {
273 if (F_marker[i]) rows_in_F[counter[i]]=i;
274 }
275 }
276 /* create mask of C nodes with value >-1, gives new id */
277 i=Paso_Util_cumsum_maskedFalse(n, mask_C, F_marker);
278 /*
279 get Prolongation :
280 */
281
282 time0=Esys_timer();
283
284 out->P=Paso_Preconditioner_AMG_getProlongation(A_p,offset_S, degree_S,S,n_C,mask_C, options->interpolation_method);
285
286 }
287
288 /*
289 construct Restriction operator as transposed of Prolongation operator:
290 */
291
292 if ( Esys_noError()) {
293 time0=Esys_timer();
294
295 out->R=Paso_Preconditioner_AMG_getRestriction(out->P);
296
297 if (SHOW_TIMING) printf("timing: level %d: Paso_SystemMatrix_getTranspose: %e\n",level,Esys_timer()-time0);
298 }
299 /*
300 construct coarse level matrix:
301 */
302 if ( Esys_noError()) {
303 time0=Esys_timer();
304
305 A_C = Paso_Preconditioner_AMG_buildInterpolationOperator(A_p, out->P, out->R);
306
307 if (SHOW_TIMING) printf("timing: level %d : construct coarse matrix: %e\n",level,Esys_timer()-time0);
308 }
309
310 /*
311 constructe courser level:
312
313 */
314 if ( Esys_noError()) {
315 out->AMG_C=Paso_Preconditioner_AMG_alloc(A_C,level+1,options);
316 }
317
318 if ( Esys_noError()) {
319 if ( out->AMG_C == NULL ) {
320 /* merge the system matrix into 1 rank when
321 it's not suitable coarsening due to the
322 total number of unknowns are too small */
323 out->A_C=A_C;
324 out->reordering = options->reordering;
325 out->refinements = options->coarse_matrix_refinements;
326 out->verbose = verbose;
327 out->options_smoother = options->smoother;
328 } else {
329 /* finally we set some helpers for the solver step */
330 out->A_C=A_C;
331 }
332 }
333 }
334 }
335 TMPMEMFREE(mask_C);
336 TMPMEMFREE(rows_in_F);
337 }
338 }
339
340 }
341 TMPMEMFREE(counter);
342 TMPMEMFREE(F_marker);
343 TMPMEMFREE(degree_S);
344 TMPMEMFREE(offset_S);
345 TMPMEMFREE(S);
346 TMPMEMFREE(degree_ST);
347 TMPMEMFREE(offset_ST);
348 TMPMEMFREE(ST);
349
350 }
351
352 if (Esys_noError()) {
353 return out;
354 } else {
355 Paso_Preconditioner_AMG_free(out);
356 return NULL;
357 }
358 }
359
360
361 void Paso_Preconditioner_AMG_solve(Paso_SystemMatrix* A, Paso_Preconditioner_AMG * amg, double * x, double * b) {
362 const dim_t n = A->mainBlock->numRows * A->mainBlock->row_block_size;
363 double time0=0;
364 const dim_t post_sweeps=amg->post_sweeps;
365 const dim_t pre_sweeps=amg->pre_sweeps;
366
367 /* presmoothing */
368 time0=Esys_timer();
369 Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, pre_sweeps, FALSE);
370
371 time0=Esys_timer()-time0;
372 if (SHOW_TIMING) printf("timing: level %d: Presmoothing: %e\n",amg->level, time0);
373 /* end of presmoothing */
374
375 if (amg->n_F < amg->n) { /* is there work on the coarse level? */
376 time0=Esys_timer();
377
378 Paso_Copy(n, amg->r, b); /* r <- b */
379 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,A,x,1.,amg->r); /*r=r-Ax*/
380 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,amg->r,0.,amg->b_C); /* b_c = R*r */
381
382 time0=Esys_timer()-time0;
383 /* coarse level solve */
384 if ( amg->AMG_C == NULL) {
385 time0=Esys_timer();
386 /* A_C is the coarsest level */
387 Paso_Preconditioner_AMG_mergeSolve(amg);
388
389 if (SHOW_TIMING) printf("timing: level %d: DIRECT SOLVER: %e\n",amg->level,Esys_timer()-time0);
390 } else {
391 Paso_Preconditioner_AMG_solve(amg->A_C, amg->AMG_C,amg->x_C,amg->b_C); /* x_C=AMG(b_C) */
392 }
393
394 time0=time0+Esys_timer();
395 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,1.,x); /* x = x + P*x_c */
396
397 /*postsmoothing*/
398
399 /*solve Ax=b with initial guess x */
400 time0=Esys_timer();
401 Paso_Preconditioner_Smoother_solve(A, amg->Smoother, x, b, post_sweeps, TRUE);
402 time0=Esys_timer()-time0;
403 if (SHOW_TIMING) printf("timing: level %d: Postsmoothing: %e\n",amg->level,time0);
404 /*end of postsmoothing*/
405 }
406
407 return;
408 }
409
410 /* theta = threshold for strong connections */
411 /* tau = threshold for diagonal dominance */
412
413 /*S_i={j \in N_i; i strongly coupled to j}
414
415 in the sense that |A_{ij}| >= theta * max_k |A_{ik}|
416 */
417
418 void Paso_Preconditioner_AMG_setStrongConnections(Paso_SystemMatrix* A,
419 dim_t *degree_S, index_t *offset_S, index_t *S,
420 const double theta, const double tau)
421 {
422
423 const dim_t my_n=A->mainBlock->numRows;
424 const dim_t overlap_n=A->row_coupleBlock->numRows;
425
426 index_t iptr, i;
427 double *threshold_p=NULL;
428
429 threshold_p = TMPMEMALLOC(2*my_n, double);
430
431 #pragma omp parallel for private(i,iptr) schedule(static)
432 for (i=0;i<my_n;++i) {
433
434 register double max_offdiagonal = 0.;
435 register double sum_row=0;
436 register double main_row=0;
437 register dim_t kdeg=0;
438 register const index_t koffset=A->mainBlock->pattern->ptr[i]+A->col_coupleBlock->pattern->ptr[i];
439
440
441 /* collect information for row i: */
442 #pragma ivdep
443 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
444 register index_t j=A->mainBlock->pattern->index[iptr];
445 register double fnorm=ABS(A->mainBlock->val[iptr]);
446 if( j != i) {
447 max_offdiagonal = MAX(max_offdiagonal,fnorm);
448 sum_row+=fnorm;
449 } else {
450 main_row=fnorm;
451 }
452
453 }
454
455 #pragma ivdep
456 for (iptr=A->col_coupleBlock->pattern->ptr[i];iptr<A->col_coupleBlock->pattern->ptr[i+1]; ++iptr) {
457 register double fnorm=ABS(A->col_coupleBlock->val[iptr]);
458
459 max_offdiagonal = MAX(max_offdiagonal,fnorm);
460 sum_row+=fnorm;
461 }
462
463 /* inspect row i: */
464 {
465 const double threshold = theta*max_offdiagonal;
466 threshold_p[2*i+1]=threshold;
467 if (tau*main_row < sum_row) { /* no diagonal dominance */
468 threshold_p[2*i]=1;
469 #pragma ivdep
470 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
471 register index_t j=A->mainBlock->pattern->index[iptr];
472 if(ABS(A->mainBlock->val[iptr])>threshold && i!=j) {
473 S[koffset+kdeg] = j;
474 kdeg++;
475 }
476 }
477 #pragma ivdep
478 for (iptr=A->col_coupleBlock->pattern->ptr[i];iptr<A->col_coupleBlock->pattern->ptr[i+1]; ++iptr) {
479 register index_t j=A->col_coupleBlock->pattern->index[iptr];
480 if(ABS(A->col_coupleBlock->val[iptr])>threshold) {
481 S[koffset+kdeg] = j + my_n;
482 kdeg++;
483 }
484 }
485 } else {
486 threshold_p[2*i]=-1;
487 }
488 }
489 offset_S[i]=koffset;
490 degree_S[i]=kdeg;
491 }
492
493 /* now we need to distribute the threshold and the diagonal dominance indicator */
494 if (A->mpi_info->size > 1) {
495
496 const index_t koffset_0=A->mainBlock->pattern->ptr[my_n]+A->col_coupleBlock->pattern->ptr[my_n]
497 -A->mainBlock->pattern->ptr[0]-A->col_coupleBlock->pattern->ptr[0];
498
499 double *remote_threshold=NULL;
500
501 Paso_Coupler* threshold_coupler=Paso_Coupler_alloc(A->row_coupler->connector ,2);
502 Paso_Coupler_startCollect(threshold_coupler,threshold_p);
503 Paso_Coupler_finishCollect(threshold_coupler);
504 remote_threshold=threshold_coupler->recv_buffer;
505
506 #pragma omp parallel for private(i,iptr) schedule(static)
507 for (i=0; i<overlap_n; i++) {
508 const double threshold = remote_threshold[2*i+1];
509 register dim_t kdeg=0;
510 register const index_t koffset=koffset_0+A->row_coupleBlock->pattern->ptr[i]+A->remote_coupleBlock->pattern->ptr[i];
511 if (remote_threshold[2*i]>0) {
512 #pragma ivdep
513 for (iptr=A->row_coupleBlock->pattern->ptr[i];iptr<A->row_coupleBlock->pattern->ptr[i+1]; ++iptr) {
514 register index_t j=A->row_coupleBlock->pattern->index[iptr];
515 if(ABS(A->row_coupleBlock->val[iptr])>threshold) {
516 S[koffset+kdeg] = j ;
517 kdeg++;
518 }
519 }
520
521 #pragma ivdep
522 for (iptr=A->remote_coupleBlock->pattern->ptr[i];iptr<A->remote_coupleBlock->pattern->ptr[i+1]; iptr++) {
523 register index_t j=A->remote_coupleBlock->pattern->index[iptr];
524 if(ABS(A->remote_coupleBlock->val[iptr])>threshold && i!=j) {
525 S[koffset+kdeg] = j + my_n;
526 kdeg++;
527 }
528 }
529 }
530 offset_S[i+my_n]=koffset;
531 degree_S[i+my_n]=kdeg;
532 }
533
534 Paso_Coupler_free(threshold_coupler);
535 }
536 TMPMEMFREE(threshold_p);
537 }
538
539 /* theta = threshold for strong connections */
540 /* tau = threshold for diagonal dominance */
541 /*S_i={j \in N_i; i strongly coupled to j}
542
543 in the sense that |A_{ij}|_F >= theta * max_k |A_{ik}|_F
544 */
545 void Paso_Preconditioner_AMG_setStrongConnections_Block(Paso_SystemMatrix* A,
546 dim_t *degree_S, index_t *offset_S, index_t *S,
547 const double theta, const double tau)
548
549 {
550 const dim_t block_size=A->block_size;
551 const dim_t my_n=A->mainBlock->numRows;
552 const dim_t overlap_n=A->row_coupleBlock->numRows;
553
554 index_t iptr, i, bi;
555 double *threshold_p=NULL;
556
557
558 threshold_p = TMPMEMALLOC(2*my_n, double);
559
560 #pragma omp parallel private(i,iptr,bi)
561 {
562
563 dim_t max_deg=0;
564 double *rtmp=NULL;
565
566 #pragma omp for schedule(static)
567 for (i=0;i<my_n;++i) max_deg=MAX(max_deg, A->mainBlock->pattern->ptr[i+1]-A->mainBlock->pattern->ptr[i]
568 +A->col_coupleBlock->pattern->ptr[i+1]-A->col_coupleBlock->pattern->ptr[i]);
569
570 rtmp=TMPMEMALLOC(max_deg, double);
571
572 #pragma omp for schedule(static)
573 for (i=0;i<my_n;++i) {
574 register double max_offdiagonal = 0.;
575 register double sum_row=0;
576 register double main_row=0;
577 register index_t rtmp_offset=-A->mainBlock->pattern->ptr[i];
578 register dim_t kdeg=0;
579 register const index_t koffset=A->mainBlock->pattern->ptr[i]+A->col_coupleBlock->pattern->ptr[i];
580
581 /* collect information for row i: */
582 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
583 register index_t j=A->mainBlock->pattern->index[iptr];
584 register double fnorm=0;
585 #pragma ivdep
586 for(bi=0;bi<block_size;++bi) {
587 register double rtmp2= A->mainBlock->val[iptr*block_size+bi];
588 fnorm+=rtmp2*rtmp2;
589 }
590 fnorm=sqrt(fnorm);
591 rtmp[iptr+rtmp_offset]=fnorm;
592
593 if( j != i) {
594 max_offdiagonal = MAX(max_offdiagonal,fnorm);
595 sum_row+=fnorm;
596 } else {
597 main_row=fnorm;
598 }
599
600 }
601
602 rtmp_offset+=A->mainBlock->pattern->ptr[i+1]-A->col_coupleBlock->pattern->ptr[i];
603 for (iptr=A->col_coupleBlock->pattern->ptr[i];iptr<A->col_coupleBlock->pattern->ptr[i+1]; ++iptr) {
604 register double fnorm=0;
605 #pragma ivdep
606 for(bi=0;bi<block_size;++bi) {
607 register double rtmp2 = A->col_coupleBlock->val[iptr*block_size+bi];
608 fnorm+=rtmp2*rtmp2;
609 }
610 fnorm=sqrt(fnorm);
611
612 rtmp[iptr+rtmp_offset]=fnorm;
613 max_offdiagonal = MAX(max_offdiagonal,fnorm);
614 sum_row+=fnorm;
615 }
616
617
618 /* inspect row i: */
619 {
620 const double threshold = theta*max_offdiagonal;
621 rtmp_offset=-A->mainBlock->pattern->ptr[i];
622
623 threshold_p[2*i+1]=threshold;
624 if (tau*main_row < sum_row) { /* no diagonal dominance */
625 threshold_p[2*i]=1;
626 #pragma ivdep
627 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
628 register index_t j=A->mainBlock->pattern->index[iptr];
629 if(rtmp[iptr+rtmp_offset] > threshold && i!=j) {
630 S[koffset+kdeg] = j;
631 kdeg++;
632 }
633 }
634 rtmp_offset+=A->mainBlock->pattern->ptr[i+1]-A->col_coupleBlock->pattern->ptr[i];
635 #pragma ivdep
636 for (iptr=A->col_coupleBlock->pattern->ptr[i];iptr<A->col_coupleBlock->pattern->ptr[i+1]; ++iptr) {
637 register index_t j=A->col_coupleBlock->pattern->index[iptr];
638 if( rtmp[iptr+rtmp_offset] >threshold) {
639 S[koffset+kdeg] = j + my_n;
640 kdeg++;
641 }
642 }
643 } else {
644 threshold_p[2*i]=-1;
645 }
646 }
647 offset_S[i]=koffset;
648 degree_S[i]=kdeg;
649 }
650 TMPMEMFREE(rtmp);
651 }
652 /* now we need to distribute the threshold and the diagonal dominance indicator */
653 if (A->mpi_info->size > 1) {
654
655 const index_t koffset_0=A->mainBlock->pattern->ptr[my_n]+A->col_coupleBlock->pattern->ptr[my_n]
656 -A->mainBlock->pattern->ptr[0]-A->col_coupleBlock->pattern->ptr[0];
657
658 double *remote_threshold=NULL;
659
660 Paso_Coupler* threshold_coupler=Paso_Coupler_alloc(A->row_coupler->connector ,2);
661 Paso_Coupler_startCollect(threshold_coupler,threshold_p);
662 Paso_Coupler_finishCollect(threshold_coupler);
663 remote_threshold=threshold_coupler->recv_buffer;
664
665 #pragma omp parallel for private(i,iptr) schedule(static)
666 for (i=0; i<overlap_n; i++) {
667
668 const double threshold2 = remote_threshold[2*i+1]*remote_threshold[2*i+1];
669 register dim_t kdeg=0;
670 register const index_t koffset=koffset_0+A->row_coupleBlock->pattern->ptr[i]+A->remote_coupleBlock->pattern->ptr[i];
671 if (remote_threshold[2*i]>0) {
672 #pragma ivdep
673 for (iptr=A->row_coupleBlock->pattern->ptr[i];iptr<A->row_coupleBlock->pattern->ptr[i+1]; ++iptr) {
674 register index_t j=A->row_coupleBlock->pattern->index[iptr];
675 register double fnorm2=0;
676 #pragma ivdepremote_threshold[2*i]
677 for(bi=0;bi<block_size;++bi) {
678 register double rtmp2 = A->row_coupleBlock->val[iptr*block_size+bi];
679 fnorm2+=rtmp2*rtmp2;
680 }
681
682 if(fnorm2 > threshold2 ) {
683 S[koffset+kdeg] = j ;
684 kdeg++;
685 }
686 }
687
688 #pragma ivdep
689 for (iptr=A->remote_coupleBlock->pattern->ptr[i];iptr<A->remote_coupleBlock->pattern->ptr[i+1]; ++iptr) {
690 register index_t j=A->remote_coupleBlock->pattern->index[iptr];
691 register double fnorm2=0;
692 #pragma ivdepremote_threshold[2*i]
693 for(bi=0;bi<block_size;++bi) {
694 register double rtmp2 = A->remote_coupleBlock->val[iptr*block_size+bi];
695 fnorm2+=rtmp2*rtmp2;
696 }
697 if(fnorm2 > threshold2 && i != j) {
698 S[koffset+kdeg] = j + my_n;
699 kdeg++;
700 }
701 }
702
703 }
704 offset_S[i+my_n]=koffset;
705 degree_S[i+my_n]=kdeg;
706 }
707 Paso_Coupler_free(threshold_coupler);
708 }
709 TMPMEMFREE(threshold_p);
710 }
711
712 void Paso_Preconditioner_AMG_transposeStrongConnections(const dim_t n, const dim_t* degree_S, const index_t* offset_S, const index_t* S,
713 const dim_t nT, dim_t* degree_ST, index_t* offset_ST,index_t* ST)
714 {
715 index_t i, j;
716 dim_t p;
717 dim_t len=0;
718 #pragma omp parallel for private(i) schedule(static)
719 for (i=0; i<nT ;++i) {
720 degree_ST[i]=0;
721 }
722 for (i=0; i<n ;++i) {
723 for (p=0; p<degree_S[i]; ++p) degree_ST[ S[offset_S[i]+p] ]++;
724 }
725 for (i=0; i<nT ;++i) {
726 offset_ST[i]=len;
727 len+=degree_ST[i];
728 degree_ST[i]=0;
729 }
730 for (i=0; i<n ;++i) {
731 for (p=0; p<degree_S[i]; ++p) {
732 j=S[offset_S[i]+p];
733 ST[offset_ST[j]+degree_ST[j]]=i;
734 degree_ST[j]++;
735 }
736 }
737 }
738
739 int compareindex(const void *a, const void *b)
740 {
741 return (*(int *)a - *(int *)b);
742 }
743
744 void Paso_Preconditioner_AMG_CIJPCoarsening(const dim_t n, const dim_t my_n, index_t*split_marker,
745 const dim_t* degree_S, const index_t* offset_S, const index_t* S,
746 const dim_t* degree_ST, const index_t* offset_ST, const index_t* ST,
747 Paso_Connector* col_connector, Paso_Distribution* col_dist)
748 {
749 dim_t i, numUndefined, iter=0;
750 index_t iptr, jptr, kptr;
751 double *random=NULL, *w=NULL, *Status=NULL;
752 index_t * ST_flag=NULL;
753
754 Paso_Coupler* w_coupler=Paso_Coupler_alloc(col_connector ,1);
755
756 w=TMPMEMALLOC(n, double);
757 Status=TMPMEMALLOC(n, double);
758 random = Paso_Distribution_createRandomVector(col_dist,1);
759 ST_flag=TMPMEMALLOC(offset_ST[n-1]+ degree_ST[n-1], index_t);
760
761 #pragma omp parallel for private(i)
762 for (i=0; i< my_n; ++i) {
763 w[i]=degree_ST[i]+random[i];
764 if (degree_ST[i] < 1) {
765 Status[i]=-100; /* F point */
766 } else {
767 Status[i]=1; /* status undefined */
768 }
769 }
770
771 #pragma omp parallel for private(i, iptr)
772 for (i=0; i< n; ++i) {
773 for( iptr =0 ; iptr < degree_ST[i]; ++iptr) {
774 ST_flag[offset_ST[i]+iptr]=1;
775 }
776 }
777
778
779 numUndefined = Paso_Distribution_numPositives(Status, col_dist, 1 );
780 /* printf(" coarsening loop start: num of undefined rows = %d \n",numUndefined); */
781 iter=0;
782 while (numUndefined > 0) {
783 Paso_Coupler_fillOverlap(n, w, w_coupler);
784
785 /* calculate the maximum value of neighbours following active strong connections:
786 w2[i]=MAX(w[k]) with k in ST[i] or S[i] and (i,k) connection is still active */
787 #pragma omp parallel for private(i, iptr)
788 for (i=0; i<my_n; ++i) {
789 if (Status[i]>0) { /* status is still undefined */
790
791 register bool_t inD=TRUE;
792 const double wi=w[i];
793
794 for( iptr =0 ; iptr < degree_S[i]; ++iptr) {
795 const index_t k=S[offset_S[i]+iptr];
796 const index_t* start_p = &ST[offset_ST[k]];
797 const index_t* where_p=(index_t*)bsearch(&i, start_p, degree_ST[k], sizeof(index_t), Paso_comparIndex);
798
799 if (ST_flag[offset_ST[k] + (index_t)(where_p-start_p)]>0) {
800 if (wi <= w[k] ) {
801 inD=FALSE;
802 break;
803 }
804 }
805
806 }
807
808 if (inD) {
809 for( iptr =0 ; iptr < degree_ST[i]; ++iptr) {
810 const index_t k=ST[offset_ST[i]+iptr];
811 if ( ST_flag[offset_ST[i]+iptr] > 0 ) {
812 if (wi <= w[k] ) {
813 inD=FALSE;
814 break;
815 }
816 }
817 }
818 }
819 if (inD) {
820 Status[i]=0.; /* is in D */
821 }
822 }
823
824 }
825
826 Paso_Coupler_fillOverlap(n, Status, w_coupler);
827
828
829 /* remove connection to D points :
830
831 for each i in D:
832 for each j in S_i:
833 w[j]--
834 ST_tag[j,i]=-1
835 for each j in ST[i]:
836 ST_tag[i,j]=-1
837 for each k in ST[j]:
838 if k in ST[i]:
839 w[j]--;
840 ST_tag[j,k]=-1
841
842 */
843 /* w is updated for local rows only */
844 {
845 #pragma omp parallel for private(i, jptr)
846 for (i=0; i< my_n; ++i) {
847
848 for (jptr=0; jptr<degree_ST[i]; ++jptr) {
849 const index_t j=ST[offset_ST[i]+jptr];
850 if ( (Status[j] == 0.) && (ST_flag[offset_ST[i]+jptr]>0) ) {
851 w[i]--;
852 ST_flag[offset_ST[i]+jptr]=-1;
853 }
854 }
855
856 }
857 #pragma omp parallel for private(i, jptr)
858 for (i=my_n; i< n; ++i) {
859 for (jptr=0; jptr<degree_ST[i]; ++jptr) {
860 const index_t j = ST[offset_ST[i]+jptr];
861 if ( Status[j] == 0. ) ST_flag[offset_ST[i]+jptr]=-1;
862 }
863 }
864
865
866 for (i=0; i< n; ++i) {
867 if ( Status[i] == 0. ) {
868
869 const index_t* start_p = &ST[offset_ST[i]];
870
871 for (jptr=0; jptr<degree_ST[i]; ++jptr) {
872 const index_t j=ST[offset_ST[i]+jptr];
873 ST_flag[offset_ST[i]+jptr]=-1;
874 for (kptr=0; kptr<degree_ST[j]; ++kptr) {
875 const index_t k=ST[offset_ST[j]+kptr];
876 if (NULL != bsearch(&k, start_p, degree_ST[i], sizeof(index_t), Paso_comparIndex) ) { /* k in ST[i] ? */
877 if (ST_flag[offset_ST[j]+kptr] >0) {
878 if (j< my_n ) {
879 w[j]--;
880 }
881 ST_flag[offset_ST[j]+kptr]=-1;
882 }
883 }
884 }
885 }
886 }
887 }
888 }
889
890 /* adjust status */
891 #pragma omp parallel for private(i)
892 for (i=0; i< my_n; ++i) {
893 if ( Status[i] == 0. ) {
894 Status[i] = -10; /* this is now a C point */
895 } else if (Status[i] == 1. && w[i]<1.) {
896 Status[i] = -100; /* this is now a F point */
897 }
898 }
899
900 i = numUndefined;
901 numUndefined = Paso_Distribution_numPositives(Status, col_dist, 1 );
902 if (numUndefined == i) {
903 Esys_setError(SYSTEM_ERROR, "Can NOT reduce numUndefined.");
904 return;
905 }
906
907 iter++;
908 /* printf(" coarsening loop %d: num of undefined rows = %d \n",iter, numUndefined); */
909
910 } /* end of while loop */
911
912 /* map to output :*/
913 Paso_Coupler_fillOverlap(n, Status, w_coupler);
914 #pragma omp parallel for private(i)
915 for (i=0; i< n; ++i) {
916 if (Status[i] > -50.) {
917 split_marker[i]=PASO_AMG_IN_C;
918 } else {
919 split_marker[i]=PASO_AMG_IN_F;
920 }
921 }
922 /* clean up : */
923 Paso_Coupler_free(w_coupler);
924 TMPMEMFREE(random);
925 TMPMEMFREE(w);
926 TMPMEMFREE(Status);
927 TMPMEMFREE(ST_flag);
928
929 return;
930 }
931
932 /* Merge the system matrix which is distributed on ranks into a complete
933 matrix on rank 0, then solve this matrix on rank 0 only */
934 Paso_SparseMatrix* Paso_Preconditioner_AMG_mergeSystemMatrix(Paso_SystemMatrix* A) {
935 index_t i, iptr, j, n, remote_n, global_n, len, offset, tag;
936 index_t row_block_size, col_block_size, block_size;
937 index_t size=A->mpi_info->size;
938 index_t rank=A->mpi_info->rank;
939 index_t *ptr=NULL, *idx=NULL, *ptr_global=NULL, *idx_global=NULL;
940 index_t *temp_n=NULL, *temp_len=NULL;
941 double *val=NULL;
942 Paso_Pattern *pattern=NULL;
943 Paso_SparseMatrix *out=NULL;
944 #ifdef ESYS_MPI
945 MPI_Request* mpi_requests=NULL;
946 MPI_Status* mpi_stati=NULL;
947 #else
948 int *mpi_requests=NULL, *mpi_stati=NULL;
949 #endif
950
951 if (size == 1) {
952 n = A->mainBlock->numRows;
953 ptr = TMPMEMALLOC(n, index_t);
954 #pragma omp parallel for private(i)
955 for (i=0; i<n; i++) ptr[i] = i;
956 out = Paso_SparseMatrix_getSubmatrix(A->mainBlock, n, n, ptr, ptr);
957 TMPMEMFREE(ptr);
958 return out;
959 }
960
961 n = A->mainBlock->numRows;
962 block_size = A->block_size;
963
964 /* Merge MainBlock and CoupleBlock to get a complete column entries
965 for each row allocated to current rank. Output (ptr, idx, val)
966 contains all info needed from current rank to merge a system
967 matrix */
968 Paso_SystemMatrix_mergeMainAndCouple(A, &ptr, &idx, &val);
969
970 #ifdef ESYS_MPI
971 mpi_requests=TMPMEMALLOC(size*2,MPI_Request);
972 mpi_stati=TMPMEMALLOC(size*2,MPI_Status);
973 #else
974 mpi_requests=TMPMEMALLOC(size*2,int);
975 mpi_stati=TMPMEMALLOC(size*2,int);
976 #endif
977
978 /* Now, pass all info to rank 0 and merge them into one sparse
979 matrix */
980 if (rank == 0) {
981 /* First, copy local ptr values into ptr_global */
982 global_n=Paso_SystemMatrix_getGlobalNumRows(A);
983 ptr_global = MEMALLOC(global_n+1, index_t);
984 memcpy(ptr_global, ptr, (n+1) * sizeof(index_t));
985 iptr = n+1;
986 MEMFREE(ptr);
987 temp_n = TMPMEMALLOC(size, index_t);
988 temp_len = TMPMEMALLOC(size, index_t);
989 temp_n[0] = iptr;
990
991 /* Second, receive ptr values from other ranks */
992 for (i=1; i<size; i++) {
993 remote_n = A->row_distribution->first_component[i+1] -
994 A->row_distribution->first_component[i];
995 #ifdef ESYS_MPI
996 MPI_Irecv(&(ptr_global[iptr]), remote_n, MPI_INT, i,
997 A->mpi_info->msg_tag_counter+i,
998 A->mpi_info->comm,
999 &mpi_requests[i]);
1000 #endif
1001 temp_n[i] = remote_n;
1002 iptr += remote_n;
1003 }
1004 #ifdef ESYS_MPI
1005 MPI_Waitall(size-1, &(mpi_requests[1]), mpi_stati);
1006 #endif
1007 A->mpi_info->msg_tag_counter += size;
1008
1009 /* Then, prepare to receive idx and val from other ranks */
1010 len = 0;
1011 offset = -1;
1012 for (i=0; i<size; i++) {
1013 if (temp_n[i] > 0) {
1014 offset += temp_n[i];
1015 len += ptr_global[offset];
1016 temp_len[i] = ptr_global[offset];
1017 }else
1018 temp_len[i] = 0;
1019 }
1020
1021 idx_global = MEMALLOC(len, index_t);
1022 iptr = temp_len[0];
1023 offset = n+1;
1024 for (i=1; i<size; i++) {
1025 len = temp_len[i];
1026 #ifdef ESYS_MPI
1027 MPI_Irecv(&(idx_global[iptr]), len, MPI_INT, i,
1028 A->mpi_info->msg_tag_counter+i,
1029 A->mpi_info->comm,
1030 &mpi_requests[i]);
1031 #endif
1032 remote_n = temp_n[i];
1033 for (j=0; j<remote_n; j++) {
1034 ptr_global[j+offset] = ptr_global[j+offset] + iptr;
1035 }
1036 offset += remote_n;
1037 iptr += len;
1038 }
1039 memcpy(idx_global, idx, temp_len[0] * sizeof(index_t));
1040 MEMFREE(idx);
1041 row_block_size = A->mainBlock->row_block_size;
1042 col_block_size = A->mainBlock->col_block_size;
1043 #ifdef ESYS_MPI
1044 MPI_Waitall(size-1, &(mpi_requests[1]), mpi_stati);
1045 #endif
1046 A->mpi_info->msg_tag_counter += size;
1047 TMPMEMFREE(temp_n);
1048
1049 /* Then generate the sparse matrix */
1050 pattern = Paso_Pattern_alloc(A->mainBlock->pattern->type, global_n,
1051 global_n, ptr_global, idx_global);
1052 out = Paso_SparseMatrix_alloc(A->mainBlock->type, pattern,
1053 row_block_size, col_block_size, FALSE);
1054 Paso_Pattern_free(pattern);
1055
1056 /* Finally, receive and copy the value */
1057 iptr = temp_len[0] * block_size;
1058 for (i=1; i<size; i++) {
1059 len = temp_len[i];
1060 #ifdef ESYS_MPI
1061 MPI_Irecv(&(out->val[iptr]), len * block_size, MPI_DOUBLE, i,
1062 A->mpi_info->msg_tag_counter+i,
1063 A->mpi_info->comm,
1064 &mpi_requests[i]);
1065 #endif
1066 iptr += (len * block_size);
1067 }
1068 memcpy(out->val, val, temp_len[0] * sizeof(double) * block_size);
1069 MEMFREE(val);
1070 #ifdef ESYS_MPI
1071 MPI_Waitall(size-1, &(mpi_requests[1]), mpi_stati);
1072 #endif
1073 A->mpi_info->msg_tag_counter += size;
1074 TMPMEMFREE(temp_len);
1075 } else { /* it's not rank 0 */
1076
1077 /* First, send out the local ptr */
1078 tag = A->mpi_info->msg_tag_counter+rank;
1079 #ifdef ESYS_MPI
1080 MPI_Issend(&(ptr[1]), n, MPI_INT, 0, tag, A->mpi_info->comm,
1081 &mpi_requests[0]);
1082 #endif
1083
1084 /* Next, send out the local idx */
1085 len = ptr[n];
1086 tag += size;
1087 #ifdef ESYS_MPI
1088 MPI_Issend(idx, len, MPI_INT, 0, tag, A->mpi_info->comm,
1089 &mpi_requests[1]);
1090 #endif
1091
1092 /* At last, send out the local val */
1093 len *= block_size;
1094 tag += size;
1095 #ifdef ESYS_MPI
1096 MPI_Issend(val, len, MPI_DOUBLE, 0, tag, A->mpi_info->comm,
1097 &mpi_requests[2]);
1098
1099 MPI_Waitall(3, mpi_requests, mpi_stati);
1100 #endif
1101 A->mpi_info->msg_tag_counter = tag + size - rank;
1102 MEMFREE(ptr);
1103 MEMFREE(idx);
1104 MEMFREE(val);
1105
1106 out = NULL;
1107 }
1108
1109 TMPMEMFREE(mpi_requests);
1110 TMPMEMFREE(mpi_stati);
1111 return out;
1112 }
1113
1114
1115 void Paso_Preconditioner_AMG_mergeSolve(Paso_Preconditioner_AMG * amg) {
1116 Paso_SystemMatrix *A = amg->A_C;
1117 Paso_SparseMatrix *A_D, *A_temp;
1118 double* x=NULL;
1119 double* b=NULL;
1120 index_t rank = A->mpi_info->rank;
1121 index_t size = A->mpi_info->size;
1122 index_t i, n, p, n_block;
1123 index_t *counts, *offset, *dist;
1124 #ifdef ESYS_MPI
1125 index_t count;
1126 #endif
1127 n_block = amg->n_block;
1128 A_D = Paso_Preconditioner_AMG_mergeSystemMatrix(A);
1129
1130 /* First, gather x and b into rank 0 */
1131 dist = A->pattern->input_distribution->first_component;
1132 n = Paso_SystemMatrix_getGlobalNumRows(A);
1133 b = TMPMEMALLOC(n*n_block, double);
1134 x = TMPMEMALLOC(n*n_block, double);
1135 counts = TMPMEMALLOC(size, index_t);
1136 offset = TMPMEMALLOC(size, index_t);
1137
1138 #pragma omp parallel for private(i,p)
1139 for (i=0; i<size; i++) {
1140 p = dist[i];
1141 counts[i] = (dist[i+1] - p)*n_block;
1142 offset[i] = p*n_block;
1143 }
1144 #ifdef ESYS_MPI
1145 count = counts[rank];
1146 MPI_Gatherv(amg->b_C, count, MPI_DOUBLE, b, counts, offset, MPI_DOUBLE, 0, A->mpi_info->comm);
1147 MPI_Gatherv(amg->x_C, count, MPI_DOUBLE, x, counts, offset, MPI_DOUBLE, 0, A->mpi_info->comm);
1148 #endif
1149
1150 if (rank == 0) {
1151 /* solve locally */
1152 #ifdef MKL
1153 A_temp = Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, A_D);
1154 A_temp->solver_package = PASO_MKL;
1155 Paso_SparseMatrix_free(A_D);
1156 Paso_MKL(A_temp, x, b, amg->reordering, amg->refinements, SHOW_TIMING);
1157 Paso_SparseMatrix_free(A_temp);
1158 #else
1159 #ifdef UMFPACK
1160 A_temp = Paso_SparseMatrix_unroll(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_CSC, A_D);
1161 A_temp->solver_package = PASO_UMFPACK;
1162 Paso_SparseMatrix_free(A_D);
1163 Paso_UMFPACK(A_temp, x, b, amg->refinements, SHOW_TIMING);
1164 Paso_SparseMatrix_free(A_temp);
1165 #else
1166 A_D->solver_p = Paso_Preconditioner_LocalSmoother_alloc(A_D, (amg->options_smoother == PASO_JACOBI), amg->verbose);
1167 A_D->solver_package = PASO_SMOOTHER;
1168 Paso_Preconditioner_LocalSmoother_solve(A_D, A_D->solver_p, x, b, amg->pre_sweeps+amg->post_sweeps, FALSE);
1169 Paso_SparseMatrix_free(A_D);
1170 #endif
1171 #endif
1172 }
1173
1174 #ifdef ESYS_MPI
1175 /* now we need to distribute the solution to all ranks */
1176 MPI_Scatterv(x, counts, offset, MPI_DOUBLE, amg->x_C, count, MPI_DOUBLE, 0, A->mpi_info->comm);
1177 #endif
1178
1179 TMPMEMFREE(x);
1180 TMPMEMFREE(b);
1181 TMPMEMFREE(counts);
1182 TMPMEMFREE(offset);
1183 }

  ViewVC Help
Powered by ViewVC 1.1.26