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

Revision 3120 - (hide annotations)
Mon Aug 30 10:48:00 2010 UTC (9 years, 7 months ago) by gross
File MIME type: text/plain
File size: 13656 byte(s)
```first iteration on Paso code clean up
```
 1 ksteube 1312 2 /******************************************************* 3 ksteube 1811 * 4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland 5 ksteube 1811 * 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 gross 495 14 ksteube 1811 15 jgs 150 /* PCG iterations */ 16 17 gross 700 #include "SystemMatrix.h" 18 #include "Paso.h" 19 jgs 150 #include "Solver.h" 20 woo409 757 21 jgs 150 #ifdef _OPENMP 22 #include 23 #endif 24 25 ksteube 1312 #ifdef PASO_MPI 26 #include 27 #endif 28 29 jgs 150 /* 30 * 31 * Purpose 32 * ======= 33 * 34 * PCG solves the linear system A*x = b using the 35 * preconditioned conjugate gradient method plus a smoother 36 * A has to be symmetric. 37 * 38 * Convergence test: norm( b - A*x )< TOL. 39 * For other measures, see the above reference. 40 * 41 * Arguments 42 * ========= 43 * 44 * r (input) DOUBLE PRECISION array, dimension N. 45 * On entry, residual of inital guess x 46 * 47 * x (input/output) DOUBLE PRECISION array, dimension N. 48 * On input, the initial guess. 49 * 50 * ITER (input/output) INT 51 * On input, the maximum iterations to be performed. 52 * On output, actual number of iterations performed. 53 * 54 * INFO (output) INT 55 * 56 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned. 57 * = SOLVEr_MAXITER_REACHED 58 * = SOLVER_INPUT_ERROR Illegal parameter: 59 * = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller 60 * = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller 61 * 62 * ============================================================== 63 */ 64 65 gross 1570 /* #define PASO_DYNAMIC_SCHEDULING_MVM */ 66 67 #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP 68 #define USE_DYNAMIC_SCHEDULING 69 #endif 70 71 jgs 150 err_t Paso_Solver_PCG( 72 Paso_SystemMatrix * A, 73 double * r, 74 double * x, 75 dim_t *iter, 76 gross 584 double * tolerance, 77 Paso_Performance* pp) { 78 jgs 150 79 /* Local variables */ 80 jfenwick 1981 dim_t num_iter=0,maxit,num_iter_global, len,rest, np, ipp; 81 #ifdef USE_DYNAMIC_SCHEDULING 82 dim_t chunk_size=-1; 83 #endif 84 gross 1570 register double ss,ss1; 85 dim_t i0, istart, iend; 86 jgs 150 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 87 err_t status = SOLVER_NO_ERROR; 88 ksteube 1312 dim_t n = Paso_SystemMatrix_getTotalNumRows(A); 89 jgs 150 double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ; 90 double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol; 91 phornby 1628 #ifdef PASO_MPI 92 double loc_sum[2], sum[2]; 93 #endif 94 jfenwick 1981 double norm_of_residual=0,norm_of_residual_global; 95 phornby 1628 register double d; 96 jgs 150 97 phornby 1628 /* Should not be any executable code before this ifdef */ 98 99 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 100 phornby 1628 101 /* Watch out for these declarations (see above) */ 102 char* chksz_chr; 103 dim_t n_chunks; 104 105 gross 1570 chksz_chr=getenv("PASO_CHUNK_SIZE_PCG"); 106 if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size); 107 phornby 1628 np=omp_get_max_threads(); 108 gross 1570 chunk_size=MIN(MAX(1,chunk_size),n/np); 109 n_chunks=n/chunk_size; 110 if (n_chunks*chunk_sizempi_info->comm); 219 ksteube 1312 #endif 220 jgs 150 tau_old=tau; 221 tau=sum_1; 222 /* p=v+beta*p */ 223 gross 1571 #pragma omp parallel private(i0, istart, iend, ipp,beta) 224 gross 1570 { 225 #ifdef USE_DYNAMIC_SCHEDULING 226 #pragma omp for schedule(dynamic, 1) 227 for (ipp=0; ipp < n_chunks; ++ipp) { 228 istart=chunk_size*ipp; 229 iend=MIN(istart+chunk_size,n); 230 #else 231 #pragma omp for schedule(static) 232 for (ipp=0; ipp mpi_info->comm); 288 ksteube 1312 #endif 289 jgs 150 delta=sum_2; 290 gross 1572 alpha=tau/delta; 291 jgs 150 292 if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) { 293 /* smoother */ 294 gross 1570 sum_3 = 0; 295 sum_4 = 0; 296 gross 1572 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1) 297 gross 1556 { 298 gross 1571 ss=0; 299 ss1=0; 300 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 301 #pragma omp for schedule(dynamic, 1) 302 for (ipp=0; ipp < n_chunks; ++ipp) { 303 istart=chunk_size*ipp; 304 iend=MIN(istart+chunk_size,n); 305 #else 306 #pragma omp for schedule(static) 307 for (ipp=0; ipp mpi_info->comm); 333 sum_3=sum[0]; 334 sum_4=sum[1]; 335 ksteube 1312 #endif 336 gross 1570 sum_5 = 0; 337 gross 1571 #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2) 338 gross 1556 { 339 jfenwick 1974 gamma_1= ( (ABS(sum_3)<= PASO_ZERO) ? 0 : -sum_4/sum_3) ; 340 gamma_2= PASO_ONE-gamma_1; 341 gross 1571 ss=0; 342 gross 1570 #ifdef USE_DYNAMIC_SCHEDULING 343 #pragma omp for schedule(dynamic, 1) 344 for (ipp=0; ipp < n_chunks; ++ipp) { 345 istart=chunk_size*ipp; 346 iend=MIN(istart+chunk_size,n); 347 #else 348 #pragma omp for schedule(static) 349 for (ipp=0; ipp mpi_info->comm); 373 ksteube 1312 #endif 374 gross 495 norm_of_residual=sqrt(sum_5); 375 convergeFlag = norm_of_residual <= tol; 376 gross 2484 maxIterFlag = num_iter > maxit; 377 gross 495 breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS); 378 } 379 gross 1556 } 380 /* end of iteration */ 381 num_iter_global=num_iter; 382 norm_of_residual_global=norm_of_residual; 383 if (maxIterFlag) { 384 status = SOLVER_MAXITER_REACHED; 385 } else if (breakFlag) { 386 status = SOLVER_BREAKDOWN; 387 } 388 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 389 jgs 150 TMPMEMFREE(rs); 390 TMPMEMFREE(x2); 391 TMPMEMFREE(v); 392 TMPMEMFREE(p); 393 *iter=num_iter_global; 394 *resid=norm_of_residual_global; 395 } 396 /* End of PCG */ 397 return status; 398 }

## Properties

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