/[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 3010 - (show annotations)
Tue Apr 27 05:10:46 2010 UTC (12 years, 3 months ago) by artak
File MIME type: text/plain
File size: 24602 byte(s)
Preparation for AMG on Systems without cut using frobenius norm
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 */
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,j;
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 /*
138 char filename[8];
139
140 sprintf(filename,"AMGLevel%d",level);
141
142 Paso_SparseMatrix_saveMM(A_p,filename);
143 */
144
145 /*Make sure we have block sizes 1*/
146 /*if (A_p->col_block_size>1) {
147 Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires column block size 1.");
148 return NULL;
149 }
150 if (A_p->row_block_size>1) {
151 Paso_setError(TYPE_ERROR,"Paso_Solver_getAMG: AMG requires row block size 1.");
152 return NULL;
153 }*/
154 out=MEMALLOC(1,Paso_Solver_AMG);
155 out->Smoother=MEMALLOC(1,Paso_Solver_Smoother);
156 /* identify independend set of rows/columns */
157 mis_marker=TMPMEMALLOC(n,index_t);
158 counter=TMPMEMALLOC(n,index_t);
159 if ( !( Paso_checkPtr(mis_marker) || Paso_checkPtr(counter) || Paso_checkPtr(out)) ) {
160 out->AMG_of_Coarse=NULL;
161 out->A_FF=NULL;
162 out->A_FC=NULL;
163 out->A_CF=NULL;
164 out->W_FC=NULL;
165 out->P=NULL;
166 out->R=NULL;
167 out->rows_in_F=NULL;
168 out->rows_in_C=NULL;
169 out->mask_F=NULL;
170 out->mask_C=NULL;
171 out->x_F=NULL;
172 out->b_F=NULL;
173 out->x_C=NULL;
174 out->b_C=NULL;
175 out->A=Paso_SparseMatrix_getReference(A_p);
176 out->solver=NULL;
177 out->Smoother->ID=options->smoother;
178 out->Smoother->Jacobi=NULL;
179 out->Smoother->GS=NULL;
180 /*out->GS=Paso_Solver_getGS(A_p,verbose);*/
181 out->level=level;
182 out->n=n;
183 out->n_F=n+1;
184 out->n_block=n_block;
185 out->post_sweeps=options->post_sweeps;
186 out->pre_sweeps=options->pre_sweeps;
187
188 sparsity=(A_p->len*1.)/(1.*A_p->numRows*A_p->numCols);
189
190 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);
191
192
193 if(sparsity>0.5) {
194 level=0;
195 }
196
197
198 if (level==0 || n<=options->min_coarse_matrix_size) {
199 out->coarsest_level=TRUE;
200 /*out->GS=Paso_Solver_getJacobi(A_p);*/
201
202 #ifdef MKL
203 out->AOffset1=Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1 + MATRIX_FORMAT_OFFSET1, out->A->pattern,1,1, FALSE);
204 #pragma omp parallel for private(i) schedule(static)
205 for (i=0;i<out->A->len;++i) {
206 out->AOffset1->val[i]=out->A->val[i];
207 }
208 #else
209 #ifdef UMFPACK
210 #else
211 if (options->smoother == PASO_JACOBI)
212 out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
213 else if (options->smoother == PASO_GS)
214 out->Smoother->GS=Paso_Solver_getGS(A_p,verbose);
215 #endif
216 #endif
217
218 } else {
219 out->coarsest_level=FALSE;
220
221 if (options->smoother == PASO_JACOBI)
222 out->Smoother->Jacobi=Paso_Solver_getJacobi(A_p);
223 else if (options->smoother == PASO_GS)
224 out->Smoother->GS=Paso_Solver_getGS(A_p,verbose);
225
226 /* identify independend set of rows/columns */
227 #pragma omp parallel for private(i) schedule(static)
228 for (i=0;i<n;++i) mis_marker[i]=-1;
229
230 /*mesuring coarsening time */
231 time0=Paso_timer();
232
233 if (options->coarsening_method == PASO_YAIR_SHAPIRA_COARSENING) {
234 Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);
235 }
236 else if (options->coarsening_method == PASO_RUGE_STUEBEN_COARSENING) {
237 Paso_Pattern_RS(A_p,mis_marker,options->coarsening_threshold);
238 }
239 else if (options->coarsening_method == PASO_AGGREGATION_COARSENING) {
240 Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);
241 }
242 else if (options->coarsening_method == PASO_STANDARD_COARSENING) {
243 Paso_Pattern_Standard_Block(A_p,mis_marker,options->coarsening_threshold);
244 }
245 else {
246 /*Default coarseneing*/
247 Paso_Pattern_Standard_Block(A_p,mis_marker,options->coarsening_threshold);
248 /*Paso_Pattern_Read("Standard.spl",n,mis_marker);*/
249 /*Paso_Pattern_YS(A_p,mis_marker,options->coarsening_threshold);*/
250 /*Paso_Pattern_greedy_Agg(A_p,mis_marker,options->coarsening_threshold);*/
251 /*Paso_Pattern_greedy(A_p->pattern,mis_marker);*/
252 /*Paso_Pattern_Aggregiation(A_p,mis_marker,options->coarsening_threshold);*/
253
254 }
255
256 if (timing) fprintf(stdout,"timing: Profilining for level %d:\n",level);
257
258 time0=Paso_timer()-time0;
259 if (timing) fprintf(stdout,"timing: Coarsening: %e\n",time0);
260
261 #pragma omp parallel for private(i) schedule(static)
262 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
263
264 out->n_F=Paso_Util_cumsum(n,counter);
265
266 if (out->n_F==0) {
267 out->coarsest_level=TRUE;
268 level=0;
269 if (verbose) {
270 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
271 printf("AMG coarsening does not eliminate any of the unknowns, switching to Jacobi preconditioner.\n");
272 }
273 }
274 else if (out->n_F==n) {
275 out->coarsest_level=TRUE;
276 level=0;
277 if (verbose) {
278 /*printf("AMG coarsening eliminates all unknowns, switching to Jacobi preconditioner.\n");*/
279 printf("AMG coarsening eliminates all of the unknowns, switching to Jacobi preconditioner.\n");
280
281 }
282 } else {
283
284 if (Paso_noError()) {
285
286 /*#pragma omp parallel for private(i) schedule(static)
287 for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
288 out->n_F=Paso_Util_cumsum(n,counter);
289 */
290
291 out->mask_F=MEMALLOC(n,index_t);
292 out->rows_in_F=MEMALLOC(out->n_F,index_t);
293 if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->rows_in_F) ) ) {
294 /* creates an index for F from mask */
295 #pragma omp parallel for private(i) schedule(static)
296 for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
297 #pragma omp parallel for private(i) schedule(static)
298 for (i = 0; i < n; ++i) {
299 if (mis_marker[i]) {
300 out->rows_in_F[counter[i]]=i;
301 out->mask_F[i]=counter[i];
302 } else {
303 out->mask_F[i]=-1;
304 }
305 }
306
307 }
308 }
309
310 /* if(level==1) {
311 printf("##TOTAL: %d, ELIMINATED: %d\n",n,out->n_F);
312 for (i = 0; i < n; ++i) {
313 printf("##%d %d\n",i,!mis_marker[i]);
314 }
315 }
316 */
317
318 /*check whether coarsening process actually makes sense to continue.
319 if coarse matrix at least smaller by 30% then continue, otherwise we stop.*/
320 if ((out->n_F*100/n)<30) {
321 level=1;
322 }
323
324 if ( Paso_noError()) {
325 out->n_C=n-out->n_F;
326 out->rows_in_C=MEMALLOC(out->n_C,index_t);
327 out->mask_C=MEMALLOC(n,index_t);
328 if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {
329 /* creates an index for C from mask */
330 #pragma omp parallel for private(i) schedule(static)
331 for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
332 Paso_Util_cumsum(n,counter);
333 #pragma omp parallel for private(i) schedule(static)
334 for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;
335 #pragma omp parallel for private(i) schedule(static)
336 for (i = 0; i < n; ++i) {
337 if (! mis_marker[i]) {
338 out->rows_in_C[counter[i]]=i;
339 out->mask_C[i]=counter[i];
340 } else {
341 out->mask_C[i]=-1;
342 }
343 }
344 }
345 }
346 if ( Paso_noError()) {
347 /* get A_FF block: */
348 /*
349 out->A_FF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_F,out->rows_in_F,out->mask_F);
350 out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
351 out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
352 */
353
354 /*Compute W_FC*/
355 /*initialy W_FC=A_FC*/
356 out->W_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
357
358 /*sprintf(filename,"W_FCbefore_%d",level);
359 Paso_SparseMatrix_saveMM(out->W_FC,filename);
360 */
361 /* for (i = 0; i < n; ++i) {
362 printf("##mis_marker[%d]=%d\n",i,mis_marker[i]);
363 }
364 */
365 time0=Paso_timer();
366 Paso_SparseMatrix_updateWeights(A_p,out->W_FC,mis_marker);
367 time0=Paso_timer()-time0;
368 if (timing) fprintf(stdout,"timing: updateWeights: %e\n",time0);
369
370 /*
371 sprintf(filename,"W_FCafter_%d",level);
372 Paso_SparseMatrix_saveMM(out->W_FC,filename);
373 */
374
375 /* get Prolongation and Restriction */
376 time0=Paso_timer();
377 out->P=Paso_SparseMatrix_getProlongation(out->W_FC,mis_marker);
378 time0=Paso_timer()-time0;
379 if (timing) fprintf(stdout,"timing: getProlongation: %e\n",time0);
380 /*out->P=Paso_SparseMatrix_loadMM_toCSR("P1.mtx");*/
381
382 /*
383 sprintf(filename,"P_%d",level);
384 Paso_SparseMatrix_saveMM(out->P,filename);
385 */
386
387 time0=Paso_timer();
388 out->R=Paso_SparseMatrix_getRestriction(out->P);
389 time0=Paso_timer()-time0;
390 if (timing) fprintf(stdout,"timing: getRestriction: %e\n",time0);
391 /*out->R=Paso_SparseMatrix_loadMM_toCSR("R1.mtx");*/
392
393 /*
394 sprintf(filename,"R_%d",level);
395 Paso_SparseMatrix_saveMM(out->R,filename);
396 */
397
398 }
399 if ( Paso_noError()) {
400
401 time0=Paso_timer();
402
403 Atemp=Paso_SparseMatrix_MatrixMatrix(A_p,out->P);
404
405 A_c=Paso_SparseMatrix_MatrixMatrix(out->R,Atemp);
406
407 /*A_c=Paso_SparseMatrix_loadMM_toCSR("A_C1.mtx");*/
408
409 Paso_SparseMatrix_free(Atemp);
410
411 /*A_c=Paso_Solver_getCoarseMatrix(A_p,out->R,out->P);*/
412 time0=Paso_timer()-time0;
413 if (timing) fprintf(stdout,"timing: getCoarseMatrix: %e\n",time0);
414
415
416 /*Paso_Solver_getCoarseMatrix(A_c, A_p,out->R,out->P);*/
417 /*
418 sprintf(filename,"A_C_%d",level);
419 Paso_SparseMatrix_saveMM(A_c,filename);
420 */
421
422 out->AMG_of_Coarse=Paso_Solver_getAMG(A_c,level-1,options);
423 }
424
425 /* allocate work arrays for AMG application */
426 if (Paso_noError()) {
427 /*
428 out->x_F=MEMALLOC(n_block*out->n_F,double);
429 out->b_F=MEMALLOC(n_block*out->n_F,double);
430 */
431 out->x_C=MEMALLOC(n_block*out->n_C,double);
432 out->b_C=MEMALLOC(n_block*out->n_C,double);
433
434 /*if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {*/
435 if ( ! ( Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {
436
437 /*
438 #pragma omp parallel for private(i) schedule(static)
439 for (i = 0; i < out->n_F; ++i) {
440 out->x_F[i]=0.;
441 out->b_F[i]=0.;
442 }
443 */
444
445 #pragma omp parallel for private(i,j) schedule(static)
446 for (i = 0; i < out->n_C; ++i) {
447 for(j=0;j<n_block;++j) {
448 out->x_C[i*n_block+j]=0.;
449 out->b_C[i*n_block+j]=0.;
450 }
451 }
452 }
453 }
454 Paso_SparseMatrix_free(A_c);
455 }
456 }
457 }
458 TMPMEMFREE(mis_marker);
459 TMPMEMFREE(counter);
460
461 if (Paso_noError()) {
462 if (verbose && level>0 && !out->coarsest_level) {
463 printf("AMG: level: %d: %d unknowns eliminated. %d left.\n",level, out->n_F,out->n_C);
464 }
465 return out;
466 } else {
467 Paso_Solver_AMG_free(out);
468 return NULL;
469 }
470 }
471
472 /**************************************************************/
473
474 /* apply AMG precondition b-> x
475
476 in fact it solves
477
478 [ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FC ] [ x_F ] = [b_F]
479 [ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C]
480
481 in the form
482
483 b->[b_F,b_C]
484 x_F=invA_FF*b_F
485 b_C=b_C-A_CF*x_F
486 x_C=AMG(b_C)
487 b_F=b_F-A_FC*x_C
488 x_F=invA_FF*b_F
489 x<-[x_F,x_C]
490
491 should be called within a parallel region
492 barrier synconization should be performed to make sure that the input vector available
493
494 */
495
496 void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) {
497 dim_t i,j;
498 double time0=0;
499 double *r=NULL, *x0=NULL;
500 bool_t timing=0;
501
502 dim_t post_sweeps=amg->post_sweeps;
503 dim_t pre_sweeps=amg->pre_sweeps;
504
505 #ifdef UMFPACK
506 Paso_UMFPACK_Handler * ptr=NULL;
507 #endif
508
509 r=MEMALLOC(amg->n*amg->n_block,double);
510 x0=MEMALLOC(amg->n*amg->n_block,double);
511
512 if (amg->coarsest_level) {
513
514 time0=Paso_timer();
515 /*If all unknown are eliminated then Jacobi is the best preconditioner*/
516
517 if (amg->n_F==0 || amg->n_F==amg->n) {
518 if(amg->Smoother->ID==PASO_JACOBI)
519 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
520 else if (amg->Smoother->ID==PASO_GS)
521 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
522 }
523 else {
524 #ifdef MKL
525 Paso_MKL1(amg->AOffset1,x,b,timing);
526 #else
527 #ifdef UMFPACK
528 ptr=(Paso_UMFPACK_Handler *)(amg->solver);
529 Paso_UMFPACK1(&ptr,amg->A,x,b,timing);
530 amg->solver=(void*) ptr;
531 #else
532 if(amg->Smoother->ID==PASO_JACOBI)
533 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
534 else if (amg->Smoother->ID==PASO_GS)
535 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
536 #endif
537 #endif
538 }
539
540 time0=Paso_timer()-time0;
541 if (timing) fprintf(stdout,"timing: DIRECT SOLVER: %e\n",time0);
542
543 } else {
544 /* presmoothing */
545 time0=Paso_timer();
546 if(amg->Smoother->ID==PASO_JACOBI)
547 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,b);
548 else if (amg->Smoother->ID==PASO_GS)
549 Paso_Solver_solveGS(amg->Smoother->GS,x,b);
550
551 /***********/
552 if (pre_sweeps>1) {
553 #pragma omp parallel for private(i,j) schedule(static)
554 for (i=0;i<amg->n;++i) {
555 for (j=0;j<amg->n_block;++j) {
556 r[i*amg->n_block+j]=b[i*amg->n_block+j];
557 }
558 }
559 }
560
561 while(pre_sweeps>1) {
562 #pragma omp parallel for private(i,j) schedule(static)
563 for (i=0;i<amg->n;++i) {
564 for (j=0;j<amg->n_block;++j) {
565 r[i*amg->n_block+j]+=b[i*amg->n_block+j];
566 }
567 }
568
569
570 /* Compute the residual r=r-Ax*/
571 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
572 /* Go round again*/
573
574 if(amg->Smoother->ID==PASO_JACOBI)
575 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x,r);
576 else if (amg->Smoother->ID==PASO_GS)
577 Paso_Solver_solveGS(amg->Smoother->GS,x,r);
578
579 pre_sweeps-=1;
580 }
581 /***********/
582
583 time0=Paso_timer()-time0;
584 if (timing) fprintf(stdout,"timing: Presmooting: %e\n",time0);
585 /* end of presmoothing */
586
587 time0=Paso_timer();
588 #pragma omp parallel for private(i,j) schedule(static)
589 for (i=0;i<amg->n;++i) {
590 for (j=0;j<amg->n_block;++j) {
591 r[i*amg->n_block+j]=b[i*amg->n_block+j];
592 }
593 }
594
595 /*r=b-Ax*/
596 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
597
598 /* b_c = R*r */
599 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->R,r,0.,amg->b_C);
600
601 time0=Paso_timer()-time0;
602 if (timing) fprintf(stdout,"timing: Before next level: %e\n",time0);
603
604 /* x_C=AMG(b_C) */
605 Paso_Solver_solveAMG(amg->AMG_of_Coarse,amg->x_C,amg->b_C);
606
607 time0=Paso_timer();
608
609 /* x_0 = P*x_c */
610 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,amg->P,amg->x_C,0.,x0);
611
612 /* x=x+x0 */
613 #pragma omp parallel for private(i,j) schedule(static)
614 for (i=0;i<amg->n;++i) {
615 for (j=0;j<amg->n_block;++j) {
616 x[i*amg->n_block+j]+=x0[i*amg->n_block+j];
617 }
618 }
619
620 /*postsmoothing*/
621
622 time0=Paso_timer();
623 #pragma omp parallel for private(i,j) schedule(static)
624 for (i=0;i<amg->n;++i) {
625 for (j=0;j<amg->n_block;++j) {
626 r[i*amg->n_block+j]=b[i*amg->n_block+j];
627 }
628 }
629
630 /*r=b-Ax */
631 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
632 if(amg->Smoother->ID==PASO_JACOBI)
633 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
634 else if (amg->Smoother->ID==PASO_GS)
635 Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
636
637 #pragma omp parallel for private(i,j) schedule(static)
638 for (i=0;i<amg->n;++i) {
639 for (j=0;j<amg->n_block;++j) {
640 x[i*amg->n_block+j]+=x0[i*amg->n_block+j];
641 }
642 }
643
644 /***************/
645 while(post_sweeps>1) {
646
647 #pragma omp parallel for private(i,j) schedule(static)
648 for (i=0;i<amg->n;++i) {
649 for (j=0;j<amg->n_block;++j) {
650 r[i*amg->n_block+j]=b[i*amg->n_block+j];
651 }
652 }
653
654 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
655
656 if(amg->Smoother->ID==PASO_JACOBI)
657 Paso_Solver_solveJacobi(amg->Smoother->Jacobi,x0,r);
658 else if (amg->Smoother->ID==PASO_GS)
659 Paso_Solver_solveGS(amg->Smoother->GS,x0,r);
660
661 #pragma omp parallel for private(i,j) schedule(static)
662 for (i=0;i<amg->n;++i) {
663 for (j=0;j<amg->n_block;++j) {
664 x[i*amg->n_block+j]+=x0[i*amg->n_block+j];
665 }
666 }
667 post_sweeps-=1;
668 }
669 /**************/
670
671 time0=Paso_timer()-time0;
672 if (timing) fprintf(stdout,"timing: Postsmoothing: %e\n",time0);
673
674 /*end of postsmoothing*/
675
676 }
677 MEMFREE(r);
678 MEMFREE(x0);
679
680 return;
681 }

  ViewVC Help
Powered by ViewVC 1.1.26