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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/paso/src/Solvers/PCG.c
File MIME type: text/plain
File size: 6312 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 /* $Id$ */
2
3 /* PCG iterations */
4
5 #include "SystemMatrix.h"
6 #include "Paso.h"
7 #include "Solver.h"
8 /* #include <math.h> */
9 #ifdef _OPENMP
10 #include <omp.h>
11 #endif
12
13 /*
14 *
15 * Purpose
16 * =======
17 *
18 * PCG solves the linear system A*x = b using the
19 * preconditioned conjugate gradient method plus a smoother
20 * A has to be symmetric.
21 *
22 * Convergence test: norm( b - A*x )< TOL.
23 * For other measures, see the above reference.
24 *
25 * Arguments
26 * =========
27 *
28 * r (input) DOUBLE PRECISION array, dimension N.
29 * On entry, residual of inital guess x
30 *
31 * x (input/output) DOUBLE PRECISION array, dimension N.
32 * On input, the initial guess.
33 *
34 * ITER (input/output) INT
35 * On input, the maximum iterations to be performed.
36 * On output, actual number of iterations performed.
37 *
38 * INFO (output) INT
39 *
40 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.
41 * = SOLVEr_MAXITER_REACHED
42 * = SOLVER_INPUT_ERROR Illegal parameter:
43 * = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller
44 * = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller
45 *
46 * ==============================================================
47 */
48
49 err_t Paso_Solver_PCG(
50 Paso_SystemMatrix * A,
51 double * r,
52 double * x,
53 dim_t *iter,
54 double * tolerance) {
55
56
57 /* Local variables */
58 dim_t num_iter=0,maxit,num_iter_global;
59 dim_t i0;
60 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
61 err_t status = SOLVER_NO_ERROR;
62 dim_t n = A->num_cols * A-> col_block_size;
63 double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;
64 double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol;
65 double d,norm_of_residual,norm_of_residual_global;
66
67 /* */
68 /*-----------------------------------------------------------------*/
69 /* */
70 /* Start of Calculation : */
71 /* --------------------- */
72 /* */
73 /* */
74 rs=TMPMEMALLOC(n,double);
75 p=TMPMEMALLOC(n,double);
76 v=TMPMEMALLOC(n,double);
77 x2=TMPMEMALLOC(n,double);
78
79 /* Test the input parameters. */
80
81 if (n < 0) {
82 status = SOLVER_INPUT_ERROR;
83 } else if (rs==NULL || p==NULL || v==NULL || x2==NULL) {
84 status = SOLVER_MEMORY_ERROR;
85 } else {
86 maxit = *iter;
87 tol = *resid;
88 #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \
89 private(tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,norm_of_residual,num_iter)
90 {
91 /* initialize data */
92 #pragma omp for private(i0) schedule(static)
93 for (i0=0;i0<n;i0++) {
94 rs[i0]=r[i0];
95 x2[i0]=x[i0];
96 p[i0]=0;
97 v[i0]=0;
98 }
99 num_iter=0;
100 /* start of iteration */
101 while (!(convergeFlag || maxIterFlag || breakFlag)) {
102 ++(num_iter);
103 #pragma omp barrier
104 #pragma omp master
105 {
106 sum_1 = 0;
107 sum_2 = 0;
108 sum_3 = 0;
109 sum_4 = 0;
110 sum_5 = 0;
111 }
112 /* v=prec(r) */
113 Paso_Solver_solvePreconditioner(A,v,r);
114 /* tau=v*r */
115 #pragma omp for private(i0) reduction(+:sum_1) schedule(static)
116 for (i0=0;i0<n;i0++) sum_1+=v[i0]*r[i0];
117 tau_old=tau;
118 tau=sum_1;
119 /* p=v+beta*p */
120 if (num_iter==1) {
121 #pragma omp for private(i0) schedule(static)
122 for (i0=0;i0<n;i0++) p[i0]=v[i0];
123 } else {
124 beta=tau/tau_old;
125 #pragma omp for private(i0) schedule(static)
126 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];
127 }
128 /* v=A*p */
129 Paso_SystemMatrix_MatrixVector(ONE, A, p,ZERO,v);
130 /* delta=p*v */
131 #pragma omp for private(i0) reduction(+:sum_2) schedule(static)
132 for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];
133 delta=sum_2;
134
135
136 if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
137 alpha=tau/delta;
138 /* smoother */
139 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)
140 for (i0=0;i0<n;i0++) {
141 r[i0]-=alpha*v[i0];
142 d=r[i0]-rs[i0];
143 sum_3=sum_3+d*d;
144 sum_4=sum_4+d*rs[i0];
145 }
146 gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
147 gamma_2= ONE-gamma_1;
148 #pragma omp for private(i0) reduction(+:sum_5) schedule(static)
149 for (i0=0;i0<n;i0++) {
150 x2[i0]+=alpha*p[i0];
151 x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
152 rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
153 sum_5+=rs[i0]*rs[i0];
154 }
155 norm_of_residual=sqrt(sum_5);
156 convergeFlag = norm_of_residual <= tol;
157 maxIterFlag = num_iter == maxit;
158 breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
159 }
160 }
161 /* end of iteration */
162 #pragma omp master
163 {
164 num_iter_global=num_iter;
165 norm_of_residual_global=norm_of_residual;
166 if (maxIterFlag) {
167 status = SOLVER_MAXITER_REACHED;
168 } else if (breakFlag) {
169 status = SOLVER_BREAKDOWN;
170 }
171 }
172 } /* end of parallel region */
173 TMPMEMFREE(rs);
174 TMPMEMFREE(x2);
175 TMPMEMFREE(v);
176 TMPMEMFREE(p);
177 *iter=num_iter_global;
178 *resid=norm_of_residual_global;
179 }
180 /* End of PCG */
181 return status;
182 }
183
184 /*
185 * $Log$
186 * Revision 1.2 2005/09/15 03:44:40 jgs
187 * Merge of development branch dev-02 back to main trunk on 2005-09-15
188 *
189 * Revision 1.1.2.1 2005/09/05 06:29:49 gross
190 * These files have been extracted from finley to define a stand alone libray for iterative
191 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
192 * has not been tested yet.
193 *
194 *
195 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26