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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2839 - (show annotations)
Wed Jan 13 23:34:49 2010 UTC (9 years, 6 months ago) by artak
File MIME type: text/plain
File size: 23451 byte(s)
minor compiliation problem under MKL is fixed.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "Solver.h"
27 #include "Options.h"
28 #include "PasoUtil.h"
29 #include "UMFPACK.h"
30 #include "MKL.h"
31 #include "SystemMatrix.h"
32 #include "Pattern_coupling.h"
33
34 /**************************************************************/
35
36 /* free all memory used by AMG */
37
38 void Paso_Solver_AMG_System_free(Paso_Solver_AMG_System * in) {
39 dim_t i;
40 if (in!=NULL) {
41 for (i=0;i<in->block_size;++i) {
42 Paso_Solver_AMG_free(in->amgblock[i]);
43 Paso_SparseMatrix_free(in->block[i]);
44 }
45 MEMFREE(in);
46 }
47 }
48
49
50 /* free all memory used by AMG */
51
52 void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {
53 if (in!=NULL) {
54
55 if(in->Smoother->ID==PASO_JACOBI)
56 Paso_Solver_Jacobi_free(in->Smoother->Jacobi);
57 else if (in->Smoother->ID==PASO_GS)
58 Paso_Solver_GS_free(in->Smoother->GS);
59 MEMFREE(in->Smoother);
60
61 Paso_SparseMatrix_free(in->A_FC);
62 Paso_SparseMatrix_free(in->A_FF);
63 Paso_SparseMatrix_free(in->W_FC);
64 Paso_SparseMatrix_free(in->A_CF);
65 Paso_SparseMatrix_free(in->P);
66 Paso_SparseMatrix_free(in->R);
67 Paso_SparseMatrix_free(in->A);
68 if(in->coarsest_level==TRUE) {
69 #ifdef MKL
70 Paso_MKL_free1(in->AOffset1);
71 Paso_SparseMatrix_free(in->AOffset1);
72 #else
73 #ifdef UMFPACK
74 Paso_UMFPACK1_free((Paso_UMFPACK_Handler*)(in->solver));
75 #endif
76 #endif
77 }
78 MEMFREE(in->rows_in_F);
79 MEMFREE(in->rows_in_C);
80 MEMFREE(in->mask_F);
81 MEMFREE(in->mask_C);
82 MEMFREE(in->x_F);
83 MEMFREE(in->b_F);
84 MEMFREE(in->x_C);
85 MEMFREE(in->b_C);
86 in->solver=NULL;
87 Paso_Solver_AMG_free(in->AMG_of_Coarse);
88 MEMFREE(in->b_C);
89 MEMFREE(in);
90 }
91 }
92
93 /**************************************************************/
94
95 /* constructs the block-block factorization of
96
97 [ A_FF A_FC ]
98 A_p=
99 [ A_CF A_FF ]
100
101 to
102
103 [ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ]
104 [ A_CF*invA_FF I ] [ 0 S ] [ 0 I ]
105
106
107 where S=A_FF-ACF*invA_FF*A_FC within the shape of S
108
109 then AMG is applied to S again until S becomes empty
110
111 */
112 Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,dim_t level,Paso_Options* options) {
113 Paso_Solver_AMG* out=NULL;
114 /*
115 Paso_Pattern* temp1=NULL;
116 Paso_Pattern* temp2=NULL;
117 */
118 bool_t verbose=options->verbose;
119 bool_t timing=0;
120
121 dim_t n=A_p->numRows;
122 dim_t n_block=A_p->row_block_size;
123 index_t* mis_marker=NULL;
124 index_t* counter=NULL;
125 /*index_t iPtr,*index, *where_p;*/
126 dim_t i;
127 Paso_SparseMatrix * A_c=NULL;
128 double time0=0;
129 Paso_SparseMatrix * Atemp=NULL;
130 double sparsity=0;
131
132 /*
133 double *temp,*temp_1;
134 double S;
135 index_t iptr;
136 */
137 /*char filename[8];*/
138 /*sprintf(filename,"AMGLevel%d",level);
139
140 Paso_SparseMatrix_saveMM(A_p,filename);
141 */
142
143 /*Make sure we have block sizes 1*/
144 if (A_p->col_block_size>1) {
145 Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1.");
146 return NULL;
147 }
148 if (A_p->row_block_size>1) {
149 Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires row block size 1.");
150 return NULL;
151 }
152 out=MEMALLOC(1,Paso_Solver_AMG);
153 out->Smoother=MEMALLOC(1,Paso_Solver_Smoother);
154 /* identify independend set of rows/columns */
155 mis_marker=TMPMEMALLOC(n,index_t);
156 counter=TMPMEMALLOC(n,index_t);
157 if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) {
158 out->AMG_of_Coarse=NULL;
159 out->A_FF=NULL;
160 out->A_FC=NULL;
161 out->A_CF=NULL;
162 out->W_FC=NULL;
163 out->P=NULL;
164 out->R=NULL;
165 out->rows_in_F=NULL;
166 out->rows_in_C=NULL;
167 out->mask_F=NULL;
168 out->mask_C=NULL;
169 out->x_F=NULL;
170 out->b_F=NULL;
171 out->x_C=NULL;
172 out->b_C=NULL;
173 out->A=Paso_SparseMatrix_getReference(A_p);
174 out->solver=NULL;
175 out->Smoother->ID=options->smoother;
176 out->Smoother->Jacobi=NULL;
177 out->Smoother->GS=NULL;
178 /*out->GS=Paso_Solver_getGS(A_p,verbose);*/
179 out->level=level;
180 out->n=n;
181 out->n_F=n+1;
182 out->n_block=n_block;
183 out->post_sweeps=options->post_sweeps;
184 out->pre_sweeps=options->pre_sweeps;
185
186 sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols);
187
188 if (verbose) fprintf(stdout,"Stats: Sparsity of the Coarse Matrix with %d non-zeros (%d,%d) in level %d is %.6f\n",A_p->len,A_p->numRows,A_p->numCols,level,sparsity);
189
190
191 if(sparsity>0.5) {
192 level=0;
193 }
194
195
196 if (level==0 || n<=options->min_coarse_matrix_size) {
197 out->coarsest_level=TRUE;
198 /*out->GS=Paso_Solver_getJacobi(A_p);*/
199
200 #ifdef MKL
201 out->AOffset1=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, out->A->pattern,1,1, FALSE);
202 #pragma omp parallel for private(i) schedule(static)
203 for (i=0;i<out->A->len;++i) {
204 out->AOffset1->val[i]=out->A->val[i];
205 }
206 #else
207 #ifdef UMFPACK
208 #else
209 if (options->smoother == PASO_JACOBI)
210 out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
211 else if (options->smoother == PASO_GS)
212 out->Smoother->GS=Paso_Solver_getGS(A_p,verbose);
213 #endif
214 #endif
215
216 } else {
217 out->coarsest_level=FALSE;
218
219 if (options->smoother == PASO_JACOBI)
220 out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
221 else if (options->smoother == PASO_GS)
222 out->Smoother->GS=Paso_Solver_getGS(A_p,verbose);
223
224 /* identify independend set of rows/columns */
225 #pragma omp parallel for private(i) schedule(static)
226 for (i=0;i<n;++i) mis_marker[i]=-1;
227
228 /*mesuring coarsening time */
229 time0=Paso_timer();
230
231 if (options->coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {
232 Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);
233 }
234 else if (options->coarsening_method == PASO_RUGE_STUEBEN_COARSENING) {
235 Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);
236 }
237 else if (options->coarsening_method == PASO_AGGREGATION_COARSENING) {
238 Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);
239 }
240 else if (options->coarsening_method == PASO_STANDARD_COARSENING) {
241 Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold);
242 }
243 else {
244 /*Default coarseneing*/
245 Paso_Pattern_Standard(A_p,mis_marker,options->coarsening_threshold);
246 /*Paso_Pattern_Read("RS.spl",n,mis_marker);*/
247 /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/
248 /*Paso_Pattern_greedy_Agg(A_p,mis_marker,options->coarsening_threshold);*/
249 /*Paso_Pattern_greedy(A_p->pattern,mis_marker);*/
250 /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
251
252 }
253
254 if (timing) fprintf(stdout,"timing: Profilining for level %d:\n",level);
255
256 time0=Paso_timer()-time0;
257 if (timing) fprintf(stdout,"timing: Coarsening: %e\n",time0);
258
259 #pragma omp parallel for private(i) schedule(static)
260 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
261
262 out->n_F=Paso_Util_cumsum(n,counter);
263
264 if (out->n_F==0) {
265 out->coarsest_level=TRUE;
266 level=0;
267 if (verbose) {
268 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
269 printf("AMG coarsening does not eliminate any of the unknowns, switching to Jacobi preconditioner.\n");
270 }
271 }
272 else if (out->n_F==n) {
273 out->coarsest_level=TRUE;
274 level=0;
275 if (verbose) {
276 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
277 printf("AMG coarsening eliminates all of the unknowns, switching to Jacobi preconditioner.\n");
278
279 }
280 } else {
281
282 if (Paso_noError()) {
283
284 /*#pragma omp parallel for private(i) schedule(static)
285 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
286 out->n_F=Paso_Util_cumsum(n,counter);
287 */
288
289 out->mask_F=MEMALLOC(n,index_t);
290 out->rows_in_F=MEMALLOC(out->n_F,index_t);
291 if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->rows_in_F) ) ) {
292 /* creates an index for F from mask */
293 #pragma omp parallel for private(i) schedule(static)
294 for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
295 #pragma omp parallel for private(i) schedule(static)
296 for (i = 0; i < n; ++i) {
297 if (mis_marker[i]) {
298 out->rows_in_F[counter[i]]=i;
299 out->mask_F[i]=counter[i];
300 } else {
301 out->mask_F[i]=-1;
302 }
303 }
304
305 }
306 }
307
308 /* if(level==1) {
309 printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F);
310 for (i = 0; i < n; ++i) {
311 printf("##%d %d\n",i,!mis_marker[i]);
312 }
313 }
314 */
315
316 /*check whether coarsening process actually makes sense to continue.
317 if coarse matrix at least smaller by 30% then continue, otherwise we stop.*/
318 if ((out->n_F*100/n)<30) {
319 level=1;
320 }
321
322 if ( Paso_noError()) {
323 out->n_C=n-out->n_F;
324 out->rows_in_C=MEMALLOC(out->n_C,index_t);
325 out->mask_C=MEMALLOC(n,index_t);
326 if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {
327 /* creates an index for C from mask */
328 #pragma omp parallel for private(i) schedule(static)
329 for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
330 Paso_Util_cumsum(n,counter);
331 #pragma omp parallel for private(i) schedule(static)
332 for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;
333 #pragma omp parallel for private(i) schedule(static)
334 for (i = 0; i < n; ++i) {
335 if (! mis_marker[i]) {
336 out->rows_in_C[counter[i]]=i;
337 out->mask_C[i]=counter[i];
338 } else {
339 out->mask_C[i]=-1;
340 }
341 }
342 }
343 }
344 if ( Paso_noError()) {
345 /* get A_FF block: */
346 /*
347 out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F);
348 out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
349 out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
350 */
351
352 /*Compute W_FC*/
353 /*initialy W_FC=A_FC*/
354 out->W_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
355
356 /*sprintf(filename,"W_FCbefore_%d",level);
357 Paso_SparseMatrix_saveMM(out->W_FC,filename);
358 */
359 /* for (i = 0; i < n; ++i) {
360 printf("##mis_marker[%d]=%d\n",i,mis_marker[i]);
361 }
362 */
363 time0=Paso_timer();
364 Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker);
365 time0=Paso_timer()-time0;
366 if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0);
367
368
369 /*sprintf(filename,"W_FCafter_%d",level);
370 Paso_SparseMatrix_saveMM(out->W_FC,filename);
371 */
372
373 /* get Prolongation and Restriction */
374 time0=Paso_timer();
375 out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker);
376 time0=Paso_timer()-time0;
377 if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0);
378 /*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/
379
380
381 /*sprintf(filename,"P_%d",level);
382 Paso_SparseMatrix_saveMM(out->P,filename);
383 */
384
385 time0=Paso_timer();
386 out->R=Paso_SparseMatrix_getRestriction(out->P);
387 time0=Paso_timer()-time0;
388 if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0);
389 /*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/
390
391 /*
392 sprintf(filename,"R_%d",level);
393 Paso_SparseMatrix_saveMM(out->R,filename);
394 */
395
396 }
397 if ( Paso_noError()) {
398
399 time0=Paso_timer();
400
401 Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
402
403 A_c=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
404
405 /*A_c=Paso_SparseMatrix_loadMM_toCSR("A_C1.mtx");*/
406
407 Paso_SparseMatrix_free(Atemp);
408
409 /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/
410 time0=Paso_timer()-time0;
411 if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0);
412
413
414 /*Paso_Solver_getCoarseMatrix(A_c, A_p,out->R,out->P);*/
415
416
417 /*sprintf(filename,"A_C_%d",level);
418 Paso_SparseMatrix_saveMM(A_c,filename);
419 */
420
421 out->AMG_of_Coarse=Paso_Solver_getAMG(A_c,level-1,options);
422 }
423
424 /* allocate work arrays for AMG application */
425 if (Paso_noError()) {
426 /*
427 out->x_F=MEMALLOC(n_block*out->n_F,double);
428 out->b_F=MEMALLOC(n_block*out->n_F,double);
429 */
430 out->x_C=MEMALLOC(n_block*out->n_C,double);
431 out->b_C=MEMALLOC(n_block*out->n_C,double);
432
433 /*if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {*/
434 if ( ! ( Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {
435
436 /*
437 #pragma omp parallel for private(i) schedule(static)
438 for (i = 0; i < out->n_F; ++i) {
439 out->x_F[i]=0.;
440 out->b_F[i]=0.;
441 }
442 */
443
444 #pragma omp parallel for private(i) schedule(static)
445 for (i = 0; i < out->n_C; ++i) {
446 out->x_C[i]=0.;
447 out->b_C[i]=0.;
448 }
449 }
450 }
451 Paso_SparseMatrix_free(A_c);
452 }
453 }
454 }
455 TMPMEMFREE(mis_marker);
456 TMPMEMFREE(counter);
457
458 if (Paso_noError()) {
459 if (verbose && level>0 && !out->coarsest_level) {
460 printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);
461 }
462 return out;
463 } else {
464 Paso_Solver_AMG_free(out);
465 return NULL;
466 }
467 }
468
469 /**************************************************************/
470
471 /* apply AMG precondition b-> x
472
473 in fact it solves
474
475 [ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FC ] [ x_F ] = [b_F]
476 [ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C]
477
478 in the form
479
480 b->[b_F,b_C]
481 x_F=invA_FF*b_F
482 b_C=b_C-A_CF*x_F
483 x_C=AMG(b_C)
484 b_F=b_F-A_FC*x_C
485 x_F=invA_FF*b_F
486 x<-[x_F,x_C]
487
488 should be called within a parallel region
489 barrier synconization should be performed to make sure that the input vector available
490
491 */
492
493 void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) {
494 dim_t i;
495 double time0=0;
496 double *r=NULL, *x0=NULL;
497 bool_t timing=0;
498
499 dim_t post_sweeps=amg->post_sweeps;
500 dim_t pre_sweeps=amg->pre_sweeps;
501
502 #ifdef UMFPACK
503 Paso_UMFPACK_Handler * ptr=NULL;
504 #endif
505
506 r=MEMALLOC(amg->n,double);
507 x0=MEMALLOC(amg->n,double);
508
509 if (amg->coarsest_level) {
510
511 time0=Paso_timer();
512 /*If all unknown are eliminated then Jacobi is the best preconditioner*/
513
514 if (amg->n_F==0 || amg->n_F==amg->n) {
515 if(amg->Smoother->ID==PASO_JACOBI)
516 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
517 else if (amg->Smoother->ID==PASO_GS)
518 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
519 }
520 else {
521 #ifdef MKL
522 Paso_MKL1(amg->AOffset1,x,b,timing);
523 #else
524 #ifdef UMFPACK
525 ptr=(Paso_UMFPACK_Handler *)(amg->solver);
526 Paso_UMFPACK1(&ptr,amg->A,x,b,timing);
527 amg->solver=(void*) ptr;
528 #else
529 if(amg->Smoother->ID==PASO_JACOBI)
530 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
531 else if (amg->Smoother->ID==PASO_GS)
532 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
533 #endif
534 #endif
535 }
536
537 time0=Paso_timer()-time0;
538 if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0);
539
540 } else {
541 /* presmoothing */
542 time0=Paso_timer();
543 if(amg->Smoother->ID==PASO_JACOBI)
544 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
545 else if (amg->Smoother->ID==PASO_GS)
546 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
547
548 /***********/
549 if (pre_sweeps>1) {
550 #pragma omp parallel for private(i) schedule(static)
551 for (i=0;i<amg->n;++i) r[i]=b[i];
552 }
553
554 while(pre_sweeps>1) {
555 #pragma omp parallel for private(i) schedule(static)
556 for (i=0;i<amg->n;++i) r[i]+=b[i];
557
558 /* Compute the residual r=r-Ax*/
559 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
560 /* Go round again*/
561
562 if(amg->Smoother->ID==PASO_JACOBI)
563 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,r);
564 else if (amg->Smoother->ID==PASO_GS)
565 Paso_Solver_solveGS(amg->Smoother->GS,x,r);
566
567 pre_sweeps-=1;
568 }
569 /***********/
570
571 time0=Paso_timer()-time0;
572 if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0);
573 /* end of presmoothing */
574
575 time0=Paso_timer();
576 #pragma omp parallel for private(i) schedule(static)
577 for (i=0;i<amg->n;++i) r[i]=b[i];
578
579 /*r=b-Ax*/
580 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
581
582 /* b_c = R*r */
583 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C);
584
585 time0=Paso_timer()-time0;
586 if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0);
587
588 /* x_C=AMG(b_C) */
589 Paso_Solver_solveAMG(amg->AMG_of_Coarse,amg->x_C,amg->b_C);
590
591 time0=Paso_timer();
592
593 /* x_0 = P*x_c */
594 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0);
595
596 /* x=x+x0 */
597 #pragma omp parallel for private(i) schedule(static)
598 for (i=0;i<amg->n;++i) x[i]+=x0[i];
599
600 /*postsmoothing*/
601
602 time0=Paso_timer();
603 #pragma omp parallel for private(i) schedule(static)
604 for (i=0;i<amg->n;++i) r[i]=b[i];
605
606 /*r=b-Ax */
607 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
608 if(amg->Smoother->ID==PASO_JACOBI)
609 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
610 else if (amg->Smoother->ID==PASO_GS)
611 Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
612
613 #pragma omp parallel for private(i) schedule(static)
614 for (i=0;i<amg->n;++i) {
615 x[i]+=x0[i];
616 }
617
618 /***************/
619 while(post_sweeps>1) {
620
621 #pragma omp parallel for private(i) schedule(static)
622 for (i=0;i<amg->n;++i) r[i]=b[i];
623
624 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
625
626 if(amg->Smoother->ID==PASO_JACOBI)
627 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
628 else if (amg->Smoother->ID==PASO_GS)
629 Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
630
631 #pragma omp parallel for private(i) schedule(static)
632 for (i=0;i<amg->n;++i) {
633 x[i]+=x0[i];
634 }
635 post_sweeps-=1;
636 }
637 /**************/
638
639 time0=Paso_timer()-time0;
640 if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);
641
642 /*end of postsmoothing*/
643
644 }
645 MEMFREE(r);
646 MEMFREE(x0);
647
648 return;
649 }

  ViewVC Help
Powered by ViewVC 1.1.26