source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_cg_relax_wt.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 10.4 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15
16/******************************************************************************
17 *
18 * ParAMG cycling routine
19 *
20 *****************************************************************************/
21
22#include "headers.h"
23#include "par_amg.h"
24
25/*--------------------------------------------------------------------------
26 * hypre_BoomerAMGCycle
27 *--------------------------------------------------------------------------*/
28
29int
30hypre_BoomerAMGCGRelaxWt( void *amg_vdata,
31 int level,
32 int num_cg_sweeps,
33 double *rlx_wt_ptr)
34{
35 hypre_ParAMGData *amg_data = amg_vdata;
36
37 MPI_Comm comm;
38 HYPRE_Solver *smoother;
39 /* Data Structure variables */
40
41 /* hypre_ParCSRMatrix **A_array = hypre_ParAMGDataAArray(amg_data); */
42 /* hypre_ParCSRMatrix **R_array = hypre_ParAMGDataRArray(amg_data); */
43 hypre_ParCSRMatrix *A = hypre_ParAMGDataAArray(amg_data)[level];
44 /* hypre_ParVector **F_array = hypre_ParAMGDataFArray(amg_data); */
45 /* hypre_ParVector **U_array = hypre_ParAMGDataUArray(amg_data); */
46 hypre_ParVector *Utemp;
47 hypre_ParVector *Vtemp;
48 hypre_ParVector *Ptemp;
49 hypre_ParVector *Rtemp;
50 hypre_ParVector *Ztemp;
51 hypre_ParVector *Qtemp;
52
53 int *CF_marker = hypre_ParAMGDataCFMarkerArray(amg_data)[level];
54 double *Ptemp_data;
55 double *Ztemp_data;
56
57 /* int **unknown_map_array;
58 int **point_map_array;
59 int **v_at_point_array; */
60
61
62 int *grid_relax_type;
63
64 /* Local variables */
65
66 int Solve_err_flag;
67 int i, j, jj;
68 int num_sweeps;
69 int relax_type;
70 int local_size;
71 int old_size;
72 int my_id = 0;
73 int smooth_type;
74 int smooth_num_levels;
75 int smooth_option = 0;
76
77 double alpha;
78 double beta;
79 double gamma = 1.0;
80 double gammaold;
81
82 double *tridiag;
83 double *trioffd;
84 double alphinv, row_sum = 0;
85 double max_row_sum = 0;
86 double rlx_wt = 0;
87 double rlx_wt_old = 0;
88 double lambda_max, lambda_max_old;
89 /* double lambda_min, lambda_min_old; */
90
91#if 0
92 double *D_mat;
93 double *S_vec;
94#endif
95
96 /* Acquire data and allocate storage */
97
98 tridiag = hypre_CTAlloc(double, num_cg_sweeps+1);
99 trioffd = hypre_CTAlloc(double, num_cg_sweeps+1);
100 for (i=0; i < num_cg_sweeps+1; i++)
101 {
102 tridiag[i] = 0;
103 trioffd[i] = 0;
104 }
105
106 Vtemp = hypre_ParAMGDataVtemp(amg_data);
107
108 Rtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
109 hypre_ParCSRMatrixGlobalNumRows(A),
110 hypre_ParCSRMatrixRowStarts(A));
111 hypre_ParVectorInitialize(Rtemp);
112 hypre_ParVectorSetPartitioningOwner(Rtemp,0);
113
114 Ptemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
115 hypre_ParCSRMatrixGlobalNumRows(A),
116 hypre_ParCSRMatrixRowStarts(A));
117 hypre_ParVectorInitialize(Ptemp);
118 hypre_ParVectorSetPartitioningOwner(Ptemp,0);
119
120 Ztemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
121 hypre_ParCSRMatrixGlobalNumRows(A),
122 hypre_ParCSRMatrixRowStarts(A));
123 hypre_ParVectorInitialize(Ztemp);
124 hypre_ParVectorSetPartitioningOwner(Ztemp,0);
125
126 Qtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
127 hypre_ParCSRMatrixGlobalNumRows(A),
128 hypre_ParCSRMatrixRowStarts(A));
129 hypre_ParVectorInitialize(Qtemp);
130 hypre_ParVectorSetPartitioningOwner(Qtemp,0);
131
132
133 grid_relax_type = hypre_ParAMGDataGridRelaxType(amg_data);
134 smooth_type = hypre_ParAMGDataSmoothType(amg_data);
135 smooth_num_levels = hypre_ParAMGDataSmoothNumLevels(amg_data);
136
137 /* Initialize */
138
139 Solve_err_flag = 0;
140
141 comm = hypre_ParCSRMatrixComm(A);
142 MPI_Comm_rank(comm,&my_id);
143
144 if (smooth_num_levels > level)
145 {
146 smoother = hypre_ParAMGDataSmoother(amg_data);
147 smooth_option = smooth_type;
148 if (smooth_type > 6 && smooth_type < 10)
149 {
150 Utemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
151 hypre_ParCSRMatrixGlobalNumRows(A),
152 hypre_ParCSRMatrixRowStarts(A));
153 hypre_ParVectorOwnsPartitioning(Utemp) = 0;
154 hypre_ParVectorInitialize(Utemp);
155 }
156 }
157
158 /*---------------------------------------------------------------------
159 * Main loop of cycling
160 *--------------------------------------------------------------------*/
161
162 relax_type = grid_relax_type[1];
163 num_sweeps = 1;
164
165 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
166 old_size
167 = hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp));
168 hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp)) =
169 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
170 Ptemp_data = hypre_VectorData(hypre_ParVectorLocalVector(Ptemp));
171 Ztemp_data = hypre_VectorData(hypre_ParVectorLocalVector(Ztemp));
172 /* if (level == 0)
173 hypre_ParVectorCopy(hypre_ParAMGDataFArray(amg_data)[0],Rtemp);
174 else
175 {
176 hypre_ParVectorCopy(F_array[level-1],Vtemp);
177 alpha = -1.0;
178 beta = 1.0;
179 hypre_ParCSRMatrixMatvec(alpha, A_array[level-1], U_array[level-1],
180 beta, Vtemp);
181 alpha = 1.0;
182 beta = 0.0;
183
184 hypre_ParCSRMatrixMatvecT(alpha,R_array[level-1],Vtemp,
185 beta,F_array[level]);
186 hypre_ParVectorCopy(F_array[level],Rtemp);
187 } */
188
189 hypre_ParVectorSetRandomValues(Rtemp,5128);
190
191 /*------------------------------------------------------------------
192 * Do the relaxation num_sweeps times
193 *-----------------------------------------------------------------*/
194
195 for (jj = 0; jj < num_cg_sweeps; jj++)
196 {
197 hypre_ParVectorSetConstantValues(Ztemp, 0.0);
198 for (j = 0; j < num_sweeps; j++)
199 {
200 {
201 Solve_err_flag = hypre_BoomerAMGRelax(A,
202 Rtemp,
203 CF_marker,
204 relax_type,
205 0,
206 1.0,
207 1.0, NULL,
208 Ztemp,
209 Vtemp, Qtemp);
210 }
211
212 if (Solve_err_flag != 0)
213 return(Solve_err_flag);
214 }
215 gammaold = gamma;
216 gamma = hypre_ParVectorInnerProd(Rtemp,Ztemp);
217 if (jj == 0)
218 {
219 hypre_ParVectorCopy(Ztemp,Ptemp);
220 beta = 1.0;
221 }
222 else
223 {
224 beta = gamma/gammaold;
225 for (i=0; i < local_size; i++)
226 Ptemp_data[i] = Ztemp_data[i] + beta*Ptemp_data[i];
227 }
228 hypre_ParCSRMatrixMatvec(1.0,A,Ptemp,0.0,Vtemp);
229 alpha = gamma /hypre_ParVectorInnerProd(Ptemp,Vtemp);
230 alphinv = 1.0/alpha;
231 tridiag[jj+1] = alphinv;
232 tridiag[jj] *= beta;
233 tridiag[jj] += alphinv;
234 trioffd[jj] *= sqrt(beta);
235 trioffd[jj+1] = -alphinv;
236 row_sum = fabs(tridiag[jj]) + fabs(trioffd[jj]);
237 if (row_sum > max_row_sum) max_row_sum = row_sum;
238 if (jj > 0)
239 {
240 row_sum = fabs(tridiag[jj-1]) + fabs(trioffd[jj-1])
241 + fabs(trioffd[jj]);
242 if (row_sum > max_row_sum) max_row_sum = row_sum;
243 /* lambda_min_old = lambda_min; */
244 lambda_max_old = lambda_max;
245 rlx_wt_old = rlx_wt;
246 hypre_Bisection(jj+1, tridiag, trioffd, lambda_max_old,
247 max_row_sum, 1.e-3, jj+1, &lambda_max);
248 rlx_wt = 1.0/lambda_max;
249 /* hypre_Bisection(jj+1, tridiag, trioffd, 0.0, lambda_min_old,
250 1.e-3, 1, &lambda_min);
251 rlx_wt = 2.0/(lambda_min+lambda_max); */
252 if (fabs(rlx_wt-rlx_wt_old) < 1.e-3 )
253 {
254 /* if (my_id == 0) printf (" cg sweeps : %d\n", (jj+1)); */
255 break;
256 }
257 }
258 else
259 {
260 /* lambda_min = tridiag[0]; */
261 lambda_max = tridiag[0];
262 }
263
264 hypre_ParVectorAxpy(-alpha,Vtemp,Rtemp);
265 }
266 /*if (my_id == 0)
267 printf (" lambda-min: %f lambda-max: %f\n", lambda_min, lambda_max);
268
269 rlx_wt = fabs(tridiag[0])+fabs(trioffd[1]);
270
271 for (i=1; i < num_cg_sweeps-1; i++)
272 {
273 row_sum = fabs(tridiag[i]) + fabs(trioffd[i]) + fabs(trioffd[i+1]);
274 if (row_sum > rlx_wt) rlx_wt = row_sum;
275 }
276 row_sum = fabs(tridiag[num_cg_sweeps-1]) + fabs(trioffd[num_cg_sweeps-1]);
277 if (row_sum > rlx_wt) rlx_wt = row_sum;
278
279 hypre_Bisection(num_cg_sweeps, tridiag, trioffd, 0.0, rlx_wt, 1.e-3, 1,
280 &lambda_min);
281 hypre_Bisection(num_cg_sweeps, tridiag, trioffd, 0.0, rlx_wt, 1.e-3,
282 num_cg_sweeps, &lambda_max);
283 */
284
285
286 hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp)) = old_size;
287
288 hypre_ParVectorDestroy(Ztemp);
289 hypre_ParVectorDestroy(Ptemp);
290 hypre_ParVectorDestroy(Rtemp);
291 hypre_ParVectorDestroy(Qtemp);
292 hypre_TFree(tridiag);
293 hypre_TFree(trioffd);
294
295 if (smooth_option > 6 && smooth_option < 10)
296 {
297 hypre_ParVectorDestroy(Utemp);
298 }
299
300 *rlx_wt_ptr = rlx_wt;
301
302 return(Solve_err_flag);
303}
304
305/*--------------------------------------------------------------------------
306 * hypre_Bisection
307 *--------------------------------------------------------------------------*/
308
309int
310hypre_Bisection(int n, double *diag, double *offd,
311 double y, double z,
312 double tol, int k, double *ev_ptr)
313{
314 double x;
315 double eigen_value;
316 int ierr = 0;
317 int sign_change = 0;
318 int i;
319 double p0, p1, p2;
320
321 while (fabs(y-z) > tol*(fabs(y) + fabs(z)))
322 {
323 x = (y+z)/2;
324
325 sign_change = 0;
326 p0 = 1;
327 p1 = diag[0] - x;
328 if (p0*p1 <= 0) sign_change++;
329 for (i=1; i < n; i++)
330 {
331 p2 = (diag[i] - x)*p1 - offd[i]*offd[i]*p0;
332 p0 = p1;
333 p1 = p2;
334 if (p0*p1 <= 0) sign_change++;
335 }
336
337 if (sign_change >= k)
338 z = x;
339 else
340 y = x;
341 }
342
343 eigen_value = (y+z)/2;
344 *ev_ptr = eigen_value;
345
346 return ierr;
347}
Note: See TracBrowser for help on using the repository browser.