# Contents of /branches/doubleplusgood/paso/src/PCG.cpp

Revision 4280 - (show annotations)
Wed Mar 6 06:45:32 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 13741 byte(s)
```some memory management replacement
```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2013 by University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development since 2012 by School of Earth Sciences 13 * 14 *****************************************************************************/ 15 16 17 /* PCG iterations */ 18 19 #include "SystemMatrix.h" 20 #include "Paso.h" 21 #include "Solver.h" 22 23 #ifdef _OPENMP 24 #include 25 #endif 26 27 #ifdef ESYS_MPI 28 #include 29 #endif 30 31 /* 32 * 33 * Purpose 34 * ======= 35 * 36 * PCG solves the linear system A*x = b using the 37 * preconditioned conjugate gradient method plus a smoother. 38 * A has to be symmetric. 39 * 40 * Convergence test: norm( b - A*x )< TOL. 41 * For other measures, see the above reference. 42 * 43 * Arguments 44 * ========= 45 * 46 * r (input) DOUBLE PRECISION array, dimension N. 47 * On entry, residual of initial guess x. 48 * 49 * x (input/output) DOUBLE PRECISION array, dimension N. 50 * On input, the initial guess. 51 * 52 * ITER (input/output) INT 53 * On input, the maximum iterations to be performed. 54 * On output, actual number of iterations performed. 55 * 56 * INFO (output) INT 57 * 58 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned. 59 * = SOLVER_MAXITER_REACHED 60 * = SOLVER_INPUT_ERROR Illegal parameter: 61 * = SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller 62 * = SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller 63 * 64 * ============================================================== 65 */ 66 67 /* #define PASO_DYNAMIC_SCHEDULING_MVM */ 68 69 #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP 70 #define USE_DYNAMIC_SCHEDULING 71 #endif 72 73 err_t Paso_Solver_PCG( 74 Paso_SystemMatrix * A, 75 double * r, 76 double * x, 77 dim_t *iter, 78 double * tolerance, 79 Paso_Performance* pp) { 80 81 /* Local variables */ 82 dim_t num_iter=0,maxit,num_iter_global, len,rest, np, ipp; 83 #ifdef USE_DYNAMIC_SCHEDULING 84 dim_t chunk_size=-1; 85 #endif 86 register double ss,ss1; 87 dim_t i0, istart, iend; 88 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 89 err_t status = SOLVER_NO_ERROR; 90 dim_t n = Paso_SystemMatrix_getTotalNumRows(A); 91 double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ; 92 double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol; 93 #ifdef ESYS_MPI 94 double loc_sum[2], sum[2]; 95 #endif 96 double norm_of_residual=0,norm_of_residual_global; 97 register double d; 98 99 /* There should not be any executable code before this ifdef */ 100 101 #ifdef USE_DYNAMIC_SCHEDULING 102 103 /* Watch out for these declarations (see above) */ 104 char* chksz_chr; 105 dim_t n_chunks; 106 107 chksz_chr=getenv("PASO_CHUNK_SIZE_PCG"); 108 if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size); 109 np=omp_get_max_threads(); 110 chunk_size=MIN(MAX(1,chunk_size),n/np); 111 n_chunks=n/chunk_size; 112 if (n_chunks*chunk_sizempi_info->comm); 221 #endif 222 tau_old=tau; 223 tau=sum_1; 224 /* p=v+beta*p */ 225 #pragma omp parallel private(i0, istart, iend, ipp,beta) 226 { 227 #ifdef USE_DYNAMIC_SCHEDULING 228 #pragma omp for schedule(dynamic, 1) 229 for (ipp=0; ipp < n_chunks; ++ipp) { 230 istart=chunk_size*ipp; 231 iend=MIN(istart+chunk_size,n); 232 #else 233 #pragma omp for schedule(static) 234 for (ipp=0; ipp mpi_info->comm); 290 #endif 291 delta=sum_2; 292 alpha=tau/delta; 293 294 if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) { 295 /* smoother */ 296 sum_3 = 0; 297 sum_4 = 0; 298 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1) 299 { 300 ss=0; 301 ss1=0; 302 #ifdef USE_DYNAMIC_SCHEDULING 303 #pragma omp for schedule(dynamic, 1) 304 for (ipp=0; ipp < n_chunks; ++ipp) { 305 istart=chunk_size*ipp; 306 iend=MIN(istart+chunk_size,n); 307 #else 308 #pragma omp for schedule(static) 309 for (ipp=0; ipp mpi_info->comm); 335 sum_3=sum[0]; 336 sum_4=sum[1]; 337 #endif 338 sum_5 = 0; 339 #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2) 340 { 341 gamma_1= ( (ABS(sum_3)<= PASO_ZERO) ? 0 : -sum_4/sum_3) ; 342 gamma_2= PASO_ONE-gamma_1; 343 ss=0; 344 #ifdef USE_DYNAMIC_SCHEDULING 345 #pragma omp for schedule(dynamic, 1) 346 for (ipp=0; ipp < n_chunks; ++ipp) { 347 istart=chunk_size*ipp; 348 iend=MIN(istart+chunk_size,n); 349 #else 350 #pragma omp for schedule(static) 351 for (ipp=0; ipp mpi_info->comm); 375 #endif 376 norm_of_residual=sqrt(sum_5); 377 convergeFlag = norm_of_residual <= tol; 378 maxIterFlag = num_iter > maxit; 379 breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS); 380 } 381 } 382 /* end of iteration */ 383 num_iter_global=num_iter; 384 norm_of_residual_global=norm_of_residual; 385 if (maxIterFlag) { 386 status = SOLVER_MAXITER_REACHED; 387 } else if (breakFlag) { 388 status = SOLVER_BREAKDOWN; 389 } 390 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 391 delete[] rs; 392 delete[] x2; 393 delete[] v; 394 delete[] p; 395 *iter=num_iter_global; 396 *resid=norm_of_residual_global; 397 } 398 /* End of PCG */ 399 return status; 400 } 401

## Properties

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