86 |
/* Local variables */ |
/* Local variables */ |
87 |
double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL, *buf1=NULL, *buf0=NULL; |
double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL, *buf1=NULL, *buf0=NULL; |
88 |
double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global; |
double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global; |
89 |
double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1; |
double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1, loc_sum[2], sum[2]; |
90 |
dim_t num_iter=0,maxit,num_iter_global; |
dim_t num_iter=0,maxit,num_iter_global; |
91 |
dim_t i0; |
dim_t i0; |
92 |
bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; |
bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; |
155 |
#pragma omp barrier |
#pragma omp barrier |
156 |
#pragma omp for private(i0) reduction(+:sum_1) schedule(static) |
#pragma omp for private(i0) reduction(+:sum_1) schedule(static) |
157 |
for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0]; |
for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0]; |
158 |
|
#ifdef PASO_MPI |
159 |
|
#pragma omp master |
160 |
|
{ |
161 |
|
loc_sum[0] = sum_1; |
162 |
|
MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
163 |
|
} |
164 |
|
#endif |
165 |
rho = sum_1; |
rho = sum_1; |
166 |
|
|
167 |
if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) { |
if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) { |
183 |
|
|
184 |
#pragma omp for private(i0) reduction(+:sum_2) schedule(static) |
#pragma omp for private(i0) reduction(+:sum_2) schedule(static) |
185 |
for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0]; |
for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0]; |
186 |
|
#ifdef PASO_MPI |
187 |
|
#pragma omp master |
188 |
|
{ |
189 |
|
loc_sum[0] = sum_2; |
190 |
|
MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
191 |
|
} |
192 |
|
#endif |
193 |
if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) { |
if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) { |
194 |
alpha = rho / sum_2; |
alpha = rho / sum_2; |
195 |
|
|
199 |
s[i0] = r[i0]; |
s[i0] = r[i0]; |
200 |
sum_3 += s[i0] * s[i0]; |
sum_3 += s[i0] * s[i0]; |
201 |
} |
} |
202 |
|
#ifdef PASO_MPI |
203 |
|
#pragma omp master |
204 |
|
{ |
205 |
|
loc_sum[0] = sum_3; |
206 |
|
MPI_Allreduce(loc_sum, &sum_3, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
207 |
|
} |
208 |
|
#endif |
209 |
norm_of_residual = sqrt(sum_3); |
norm_of_residual = sqrt(sum_3); |
210 |
|
|
211 |
/* Early check for tolerance. */ |
/* Early check for tolerance. */ |
224 |
omegaNumtr +=t[i0] * s[i0]; |
omegaNumtr +=t[i0] * s[i0]; |
225 |
omegaDenumtr += t[i0] * t[i0]; |
omegaDenumtr += t[i0] * t[i0]; |
226 |
} |
} |
227 |
|
#ifdef PASO_MPI |
228 |
|
#pragma omp master |
229 |
|
{ |
230 |
|
loc_sum[0] = omegaNumtr; |
231 |
|
loc_sum[1] = omegaDenumtr; |
232 |
|
MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
233 |
|
omegaNumtr=sum[0]; |
234 |
|
omegaDenumtr=sum[1]; |
235 |
|
} |
236 |
|
#endif |
237 |
if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) { |
if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) { |
238 |
omega = omegaNumtr / omegaDenumtr; |
omega = omegaNumtr / omegaDenumtr; |
239 |
|
|
243 |
r[i0] = s[i0]-omega * t[i0]; |
r[i0] = s[i0]-omega * t[i0]; |
244 |
sum_4 += r[i0] * r[i0]; |
sum_4 += r[i0] * r[i0]; |
245 |
} |
} |
246 |
|
#ifdef PASO_MPI |
247 |
|
#pragma omp master |
248 |
|
{ |
249 |
|
loc_sum[0] = sum_4; |
250 |
|
MPI_Allreduce(loc_sum, &sum_4, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
251 |
|
} |
252 |
|
#endif |
253 |
norm_of_residual = sqrt(sum_4); |
norm_of_residual = sqrt(sum_4); |
254 |
convergeFlag = norm_of_residual <= tol; |
convergeFlag = norm_of_residual <= tol; |
255 |
maxIterFlag = num_iter == maxit; |
maxIterFlag = num_iter == maxit; |