1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * ParAMG cycling routine
11 *
12 *****************************************************************************/
13
14 #include "_hypre_parcsr_ls.h"
15 #include "par_amg.h"
16
17 /*--------------------------------------------------------------------------
18 * hypre_BoomerAMGCycle
19 *--------------------------------------------------------------------------*/
20
21 HYPRE_Int
hypre_BoomerAMGCGRelaxWt(void * amg_vdata,HYPRE_Int level,HYPRE_Int num_cg_sweeps,HYPRE_Real * rlx_wt_ptr)22 hypre_BoomerAMGCGRelaxWt( void *amg_vdata,
23 HYPRE_Int level,
24 HYPRE_Int num_cg_sweeps,
25 HYPRE_Real *rlx_wt_ptr)
26 {
27 hypre_ParAMGData *amg_data = (hypre_ParAMGData*) amg_vdata;
28
29 MPI_Comm comm;
30 HYPRE_Solver *smoother;
31 /* Data Structure variables */
32
33 /* hypre_ParCSRMatrix **A_array = hypre_ParAMGDataAArray(amg_data); */
34 /* hypre_ParCSRMatrix **R_array = hypre_ParAMGDataRArray(amg_data); */
35 hypre_ParCSRMatrix *A = hypre_ParAMGDataAArray(amg_data)[level];
36 /* hypre_ParVector **F_array = hypre_ParAMGDataFArray(amg_data); */
37 /* hypre_ParVector **U_array = hypre_ParAMGDataUArray(amg_data); */
38 hypre_ParVector *Utemp;
39 hypre_ParVector *Vtemp;
40 hypre_ParVector *Ptemp;
41 hypre_ParVector *Rtemp;
42 hypre_ParVector *Ztemp;
43 hypre_ParVector *Qtemp = NULL;
44
45 HYPRE_Int *CF_marker = hypre_IntArrayData(hypre_ParAMGDataCFMarkerArray(amg_data)[level]);
46 HYPRE_Real *Ptemp_data;
47 HYPRE_Real *Ztemp_data;
48
49 /* HYPRE_Int **unknown_map_array;
50 HYPRE_Int **point_map_array;
51 HYPRE_Int **v_at_point_array; */
52
53
54 HYPRE_Int *grid_relax_type;
55
56 /* Local variables */
57 HYPRE_Int Solve_err_flag;
58 HYPRE_Int i, j, jj;
59 HYPRE_Int num_sweeps;
60 HYPRE_Int relax_type;
61 HYPRE_Int local_size;
62 HYPRE_Int old_size;
63 HYPRE_Int my_id = 0;
64 HYPRE_Int smooth_type;
65 HYPRE_Int smooth_num_levels;
66 HYPRE_Int smooth_option = 0;
67 HYPRE_Int needQ = 0;
68
69 hypre_Vector *l1_norms = NULL;
70
71 HYPRE_Real alpha;
72 HYPRE_Real beta;
73 HYPRE_Real gamma = 1.0;
74 HYPRE_Real gammaold;
75
76 HYPRE_Real *tridiag;
77 HYPRE_Real *trioffd;
78 HYPRE_Real alphinv, row_sum = 0;
79 HYPRE_Real max_row_sum = 0;
80 HYPRE_Real rlx_wt = 0;
81 HYPRE_Real rlx_wt_old = 0;
82 HYPRE_Real lambda_max, lambda_max_old;
83 /* HYPRE_Real lambda_min, lambda_min_old; */
84
85 #if 0
86 HYPRE_Real *D_mat;
87 HYPRE_Real *S_vec;
88 #endif
89
90 #if !defined(HYPRE_USING_CUDA) && !defined(HYPRE_USING_HIP)
91 HYPRE_Int num_threads = hypre_NumThreads();
92 #endif
93
94 /* Acquire data and allocate storage */
95
96 tridiag = hypre_CTAlloc(HYPRE_Real, num_cg_sweeps+1, HYPRE_MEMORY_HOST);
97 trioffd = hypre_CTAlloc(HYPRE_Real, num_cg_sweeps+1, HYPRE_MEMORY_HOST);
98 for (i=0; i < num_cg_sweeps+1; i++)
99 {
100 tridiag[i] = 0;
101 trioffd[i] = 0;
102 }
103
104 Vtemp = hypre_ParAMGDataVtemp(amg_data);
105
106 Rtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
107 hypre_ParCSRMatrixGlobalNumRows(A),
108 hypre_ParCSRMatrixRowStarts(A));
109 hypre_ParVectorInitialize(Rtemp);
110
111 Ptemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
112 hypre_ParCSRMatrixGlobalNumRows(A),
113 hypre_ParCSRMatrixRowStarts(A));
114 hypre_ParVectorInitialize(Ptemp);
115
116 Ztemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
117 hypre_ParCSRMatrixGlobalNumRows(A),
118 hypre_ParCSRMatrixRowStarts(A));
119 hypre_ParVectorInitialize(Ztemp);
120
121 if (hypre_ParAMGDataL1Norms(amg_data) != NULL)
122 l1_norms = hypre_ParAMGDataL1Norms(amg_data)[level];
123
124 #if !defined(HYPRE_USING_CUDA) && !defined(HYPRE_USING_HIP)
125 if (num_threads > 1)
126 #endif
127 {
128 needQ = 1;
129 }
130
131 grid_relax_type = hypre_ParAMGDataGridRelaxType(amg_data);
132 smooth_type = hypre_ParAMGDataSmoothType(amg_data);
133 smooth_num_levels = hypre_ParAMGDataSmoothNumLevels(amg_data);
134
135 /* Initialize */
136
137 Solve_err_flag = 0;
138
139 comm = hypre_ParCSRMatrixComm(A);
140 hypre_MPI_Comm_rank(comm,&my_id);
141
142 if (smooth_num_levels > level)
143 {
144 smoother = hypre_ParAMGDataSmoother(amg_data);
145 smooth_option = smooth_type;
146 if (smooth_type > 6 && smooth_type < 10)
147 {
148 Utemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
149 hypre_ParCSRMatrixGlobalNumRows(A),
150 hypre_ParCSRMatrixRowStarts(A));
151 hypre_ParVectorInitialize(Utemp);
152 }
153 }
154
155 /*---------------------------------------------------------------------
156 * Main loop of cycling
157 *--------------------------------------------------------------------*/
158
159 relax_type = grid_relax_type[1];
160 num_sweeps = 1;
161
162 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
163 old_size = hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp));
164 hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp)) =
165 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
166 Ptemp_data = hypre_VectorData(hypre_ParVectorLocalVector(Ptemp));
167 Ztemp_data = hypre_VectorData(hypre_ParVectorLocalVector(Ztemp));
168 /* if (level == 0)
169 hypre_ParVectorCopy(hypre_ParAMGDataFArray(amg_data)[0],Rtemp);
170 else
171 {
172 hypre_ParVectorCopy(F_array[level-1],Vtemp);
173 alpha = -1.0;
174 beta = 1.0;
175 hypre_ParCSRMatrixMatvec(alpha, A_array[level-1], U_array[level-1],
176 beta, Vtemp);
177 alpha = 1.0;
178 beta = 0.0;
179
180 hypre_ParCSRMatrixMatvecT(alpha,R_array[level-1],Vtemp,
181 beta,F_array[level]);
182 hypre_ParVectorCopy(F_array[level],Rtemp);
183 } */
184
185 hypre_ParVectorSetRandomValues(Rtemp, 5128);
186
187 if (needQ)
188 {
189 Qtemp = hypre_ParMultiVectorCreate(hypre_ParCSRMatrixComm(A),
190 hypre_ParCSRMatrixGlobalNumRows(A),
191 hypre_ParCSRMatrixRowStarts(A),
192 needQ);
193 hypre_ParVectorInitialize(Qtemp);
194 }
195
196 /*------------------------------------------------------------------
197 * Do the relaxation num_sweeps times
198 *-----------------------------------------------------------------*/
199
200 for (jj = 0; jj < num_cg_sweeps; jj++)
201 {
202 hypre_ParVectorSetConstantValues(Ztemp, 0.0);
203 for (j = 0; j < num_sweeps; j++)
204 {
205 if (smooth_option > 6)
206 {
207
208 hypre_ParVectorCopy(Rtemp,Vtemp);
209 alpha = -1.0;
210 beta = 1.0;
211 hypre_ParCSRMatrixMatvec(alpha, A,
212 Ztemp, beta, Vtemp);
213 if (smooth_option == 8)
214 HYPRE_ParCSRParaSailsSolve(smoother[level],
215 (HYPRE_ParCSRMatrix) A,
216 (HYPRE_ParVector) Vtemp,
217 (HYPRE_ParVector) Utemp);
218 else if (smooth_option == 7)
219 {
220 HYPRE_ParCSRPilutSolve(smoother[level],
221 (HYPRE_ParCSRMatrix) A,
222 (HYPRE_ParVector) Vtemp,
223 (HYPRE_ParVector) Utemp);
224 hypre_ParVectorAxpy(1.0,Utemp,Ztemp);
225 }
226 else if (smooth_option == 9)
227 {
228 HYPRE_EuclidSolve(smoother[level],
229 (HYPRE_ParCSRMatrix) A,
230 (HYPRE_ParVector) Vtemp,
231 (HYPRE_ParVector) Utemp);
232 hypre_ParVectorAxpy(1.0,Utemp,Ztemp);
233 }
234 }
235 else if (smooth_option == 6)
236 HYPRE_SchwarzSolve(smoother[level],
237 (HYPRE_ParCSRMatrix) A,
238 (HYPRE_ParVector) Rtemp,
239 (HYPRE_ParVector) Ztemp);
240 else
241 {
242 Solve_err_flag = hypre_BoomerAMGRelax(A,
243 Rtemp,
244 CF_marker,
245 relax_type,
246 0,
247 1.0,
248 1.0,
249 l1_norms ? hypre_VectorData(l1_norms) : NULL,
250 Ztemp,
251 Vtemp,
252 Qtemp);
253 }
254
255 if (Solve_err_flag != 0)
256 {
257 hypre_ParVectorDestroy(Ptemp);
258 hypre_TFree(tridiag, HYPRE_MEMORY_HOST);
259 hypre_TFree(trioffd, HYPRE_MEMORY_HOST);
260 return(Solve_err_flag);
261 }
262 }
263 gammaold = gamma;
264 gamma = hypre_ParVectorInnerProd(Rtemp,Ztemp);
265 if (jj == 0)
266 {
267 hypre_ParVectorCopy(Ztemp,Ptemp);
268 beta = 1.0;
269 }
270 else
271 {
272 beta = gamma/gammaold;
273 for (i=0; i < local_size; i++)
274 Ptemp_data[i] = Ztemp_data[i] + beta*Ptemp_data[i];
275 }
276 hypre_ParCSRMatrixMatvec(1.0,A,Ptemp,0.0,Vtemp);
277 alpha = gamma /hypre_ParVectorInnerProd(Ptemp,Vtemp);
278 alphinv = 1.0/alpha;
279 tridiag[jj+1] = alphinv;
280 tridiag[jj] *= beta;
281 tridiag[jj] += alphinv;
282 trioffd[jj] *= sqrt(beta);
283 trioffd[jj+1] = -alphinv;
284 row_sum = fabs(tridiag[jj]) + fabs(trioffd[jj]);
285 if (row_sum > max_row_sum) max_row_sum = row_sum;
286 if (jj > 0)
287 {
288 row_sum = fabs(tridiag[jj-1]) + fabs(trioffd[jj-1])
289 + fabs(trioffd[jj]);
290 if (row_sum > max_row_sum) max_row_sum = row_sum;
291 /* lambda_min_old = lambda_min; */
292 lambda_max_old = lambda_max;
293 rlx_wt_old = rlx_wt;
294 hypre_Bisection(jj+1, tridiag, trioffd, lambda_max_old,
295 max_row_sum, 1.e-3, jj+1, &lambda_max);
296 rlx_wt = 1.0/lambda_max;
297 /* hypre_Bisection(jj+1, tridiag, trioffd, 0.0, lambda_min_old,
298 1.e-3, 1, &lambda_min);
299 rlx_wt = 2.0/(lambda_min+lambda_max); */
300 if (fabs(rlx_wt-rlx_wt_old) < 1.e-3 )
301 {
302 /* if (my_id == 0) hypre_printf (" cg sweeps : %d\n", (jj+1)); */
303 break;
304 }
305 }
306 else
307 {
308 /* lambda_min = tridiag[0]; */
309 lambda_max = tridiag[0];
310 }
311
312 hypre_ParVectorAxpy(-alpha,Vtemp,Rtemp);
313 }
314 /*if (my_id == 0)
315 hypre_printf (" lambda-min: %f lambda-max: %f\n", lambda_min, lambda_max);
316
317 rlx_wt = fabs(tridiag[0])+fabs(trioffd[1]);
318
319 for (i=1; i < num_cg_sweeps-1; i++)
320 {
321 row_sum = fabs(tridiag[i]) + fabs(trioffd[i]) + fabs(trioffd[i+1]);
322 if (row_sum > rlx_wt) rlx_wt = row_sum;
323 }
324 row_sum = fabs(tridiag[num_cg_sweeps-1]) + fabs(trioffd[num_cg_sweeps-1]);
325 if (row_sum > rlx_wt) rlx_wt = row_sum;
326
327 hypre_Bisection(num_cg_sweeps, tridiag, trioffd, 0.0, rlx_wt, 1.e-3, 1,
328 &lambda_min);
329 hypre_Bisection(num_cg_sweeps, tridiag, trioffd, 0.0, rlx_wt, 1.e-3,
330 num_cg_sweeps, &lambda_max);
331 */
332
333
334 hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp)) = old_size;
335
336 hypre_ParVectorDestroy(Ztemp);
337 hypre_ParVectorDestroy(Ptemp);
338 hypre_ParVectorDestroy(Rtemp);
339
340 if (Qtemp)
341 {
342 hypre_ParVectorDestroy(Qtemp);
343 }
344
345 hypre_TFree(tridiag, HYPRE_MEMORY_HOST);
346 hypre_TFree(trioffd, HYPRE_MEMORY_HOST);
347
348 if (smooth_option > 6 && smooth_option < 10)
349 {
350 hypre_ParVectorDestroy(Utemp);
351 }
352
353 *rlx_wt_ptr = rlx_wt;
354
355 return(Solve_err_flag);
356 }
357
358 /*--------------------------------------------------------------------------
359 * hypre_Bisection
360 *--------------------------------------------------------------------------*/
361
362 HYPRE_Int
hypre_Bisection(HYPRE_Int n,HYPRE_Real * diag,HYPRE_Real * offd,HYPRE_Real y,HYPRE_Real z,HYPRE_Real tol,HYPRE_Int k,HYPRE_Real * ev_ptr)363 hypre_Bisection(HYPRE_Int n, HYPRE_Real *diag, HYPRE_Real *offd,
364 HYPRE_Real y, HYPRE_Real z,
365 HYPRE_Real tol, HYPRE_Int k, HYPRE_Real *ev_ptr)
366 {
367 HYPRE_Real x;
368 HYPRE_Real eigen_value;
369 HYPRE_Int ierr = 0;
370 HYPRE_Int sign_change = 0;
371 HYPRE_Int i;
372 HYPRE_Real p0, p1, p2;
373
374 while (fabs(y-z) > tol*(fabs(y) + fabs(z)))
375 {
376 x = (y+z)/2;
377
378 sign_change = 0;
379 p0 = 1;
380 p1 = diag[0] - x;
381 if (p0*p1 <= 0) sign_change++;
382 for (i=1; i < n; i++)
383 {
384 p2 = (diag[i] - x)*p1 - offd[i]*offd[i]*p0;
385 p0 = p1;
386 p1 = p2;
387 if (p0*p1 <= 0) sign_change++;
388 }
389
390 if (sign_change >= k)
391 z = x;
392 else
393 y = x;
394 }
395
396 eigen_value = (y+z)/2;
397 *ev_ptr = eigen_value;
398
399 return ierr;
400 }
401