# Annotation of /trunk/paso/src/PCG.c

Revision 1572 - (hide annotations)
Mon May 26 12:50:56 2008 UTC (14 years, 2 months ago) by gross
File MIME type: text/plain
File size: 13207 byte(s)
```some openmp problems in PCG fixed.
```
 1 ksteube 1312 2 jgs 150 /* \$Id\$ */ 3 4 ksteube 1312 /******************************************************* 5 * 6 * Copyright 2003-2007 by ACceSS MNRF 7 * Copyright 2007 by University of Queensland 8 * 9 * http://esscc.uq.edu.au 10 * Primary Business: Queensland, Australia 11 * Licensed under the Open Software License version 3.0 12 * http://www.opensource.org/licenses/osl-3.0.php 13 * 14 *******************************************************/ 15 gross 495 16 jgs 150 /* PCG iterations */ 17 18 gross 700 #include "SystemMatrix.h" 19 #include "Paso.h" 20 jgs 150 #include "Solver.h" 21 woo409 757 22 jgs 150 #ifdef _OPENMP 23 #include 24 #endif 25 26 ksteube 1312 #ifdef PASO_MPI 27 #include 28 #endif 29 30 jgs 150 /* 31 * 32 * Purpose 33 * ======= 34 * 35 * PCG solves the linear system A*x = b using the 36 * preconditioned conjugate gradient method plus a smoother 37 * A has to be symmetric. 38 * 39 * Convergence test: norm( b - A*x )< TOL. 40 * For other measures, see the above reference. 41 * 42 * Arguments 43 * ========= 44 * 45 * r (input) DOUBLE PRECISION array, dimension N. 46 * On entry, residual of inital guess x 47 * 48 * x (input/output) DOUBLE PRECISION array, dimension N. 49 * On input, the initial guess. 50 * 51 * ITER (input/output) INT 52 * On input, the maximum iterations to be performed. 53 * On output, actual number of iterations performed. 54 * 55 * INFO (output) INT 56 * 57 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned. 58 * = SOLVEr_MAXITER_REACHED 59 * = SOLVER_INPUT_ERROR Illegal parameter: 60 * = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller 61 * = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller 62 * 63 * ============================================================== 64 */ 65 66 gross 1570 /* #define PASO_DYNAMIC_SCHEDULING_MVM */ 67 68 #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP 69 #define USE_DYNAMIC_SCHEDULING 70 #endif 71 72 jgs 150 err_t Paso_Solver_PCG( 73 Paso_SystemMatrix * A, 74 double * r, 75 double * x, 76 dim_t *iter, 77 gross 584 double * tolerance, 78 Paso_Performance* pp) { 79 jgs 150 80 /* Local variables */ 81 gross 1570 dim_t num_iter=0,maxit,num_iter_global, chunk_size=-1, len,rest, n_chunks, np, ipp; 82 register double ss,ss1; 83 dim_t i0, istart, iend; 84 jgs 150 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 85 err_t status = SOLVER_NO_ERROR; 86 ksteube 1312 dim_t n = Paso_SystemMatrix_getTotalNumRows(A); 87 jgs 150 double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ; 88 double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol; 89 ksteube 1312 double norm_of_residual,norm_of_residual_global, loc_sum[2], sum[2]; 90 gross 495 register double r_tmp,d,rs_tmp,x2_tmp,x_tmp; 91 gross 1570 char* chksz_chr; 92 np=omp_get_max_threads(); 93 jgs 150 94 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 95 chksz_chr=getenv("PASO_CHUNK_SIZE_PCG"); 96 if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size); 97 chunk_size=MIN(MAX(1,chunk_size),n/np); 98 n_chunks=n/chunk_size; 99 if (n_chunks*chunk_sizempi_info->comm); 196 ksteube 1312 #endif 197 jgs 150 tau_old=tau; 198 tau=sum_1; 199 /* p=v+beta*p */ 200 gross 1571 #pragma omp parallel private(i0, istart, iend, ipp,beta) 201 gross 1570 { 202 #ifdef USE_DYNAMIC_SCHEDULING 203 #pragma omp for schedule(dynamic, 1) 204 for (ipp=0; ipp < n_chunks; ++ipp) { 205 istart=chunk_size*ipp; 206 iend=MIN(istart+chunk_size,n); 207 #else 208 #pragma omp for schedule(static) 209 for (ipp=0; ipp mpi_info->comm); 265 ksteube 1312 #endif 266 jgs 150 delta=sum_2; 267 gross 1572 alpha=tau/delta; 268 jgs 150 269 if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) { 270 /* smoother */ 271 gross 1570 sum_3 = 0; 272 sum_4 = 0; 273 gross 1572 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1) 274 gross 1556 { 275 gross 1571 ss=0; 276 ss1=0; 277 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 278 #pragma omp for schedule(dynamic, 1) 279 for (ipp=0; ipp < n_chunks; ++ipp) { 280 istart=chunk_size*ipp; 281 iend=MIN(istart+chunk_size,n); 282 #else 283 #pragma omp for schedule(static) 284 for (ipp=0; ipp mpi_info->comm); 310 sum_3=sum[0]; 311 sum_4=sum[1]; 312 ksteube 1312 #endif 313 gross 1570 sum_5 = 0; 314 gross 1571 #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2) 315 gross 1556 { 316 gross 1570 gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ; 317 gamma_2= ONE-gamma_1; 318 gross 1571 ss=0; 319 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 320 #pragma omp for schedule(dynamic, 1) 321 for (ipp=0; ipp < n_chunks; ++ipp) { 322 istart=chunk_size*ipp; 323 iend=MIN(istart+chunk_size,n); 324 #else 325 #pragma omp for schedule(static) 326 for (ipp=0; ipp mpi_info->comm); 350 ksteube 1312 #endif 351 gross 495 norm_of_residual=sqrt(sum_5); 352 convergeFlag = norm_of_residual <= tol; 353 maxIterFlag = num_iter == maxit; 354 breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS); 355 } 356 gross 1556 } 357 /* end of iteration */ 358 num_iter_global=num_iter; 359 norm_of_residual_global=norm_of_residual; 360 if (maxIterFlag) { 361 status = SOLVER_MAXITER_REACHED; 362 } else if (breakFlag) { 363 status = SOLVER_BREAKDOWN; 364 } 365 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 366 jgs 150 TMPMEMFREE(rs); 367 TMPMEMFREE(x2); 368 TMPMEMFREE(v); 369 TMPMEMFREE(p); 370 *iter=num_iter_global; 371 *resid=norm_of_residual_global; 372 } 373 /* End of PCG */ 374 return status; 375 }

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision