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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 1 month 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 jgs 150 /* $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