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 #include "_hypre_sstruct_ls.h"
9
10 /*--------------------------------------------------------------------------
11 * hypre_MaxwellSolve- note that there is no input operator Aee. We assume
12 * that maxwell_vdata has the exact operators. This prevents the need to
13 * to recompute Ann in the solve phase. However, we do allow the f_edge &
14 * u_edge to change per call.
15 *--------------------------------------------------------------------------*/
16
17 HYPRE_Int
hypre_MaxwellSolve(void * maxwell_vdata,hypre_SStructMatrix * A_in,hypre_SStructVector * f,hypre_SStructVector * u)18 hypre_MaxwellSolve( void *maxwell_vdata,
19 hypre_SStructMatrix *A_in,
20 hypre_SStructVector *f,
21 hypre_SStructVector *u )
22 {
23 hypre_MaxwellData *maxwell_data = (hypre_MaxwellData *) maxwell_vdata;
24
25 hypre_ParVector *f_edge;
26 hypre_ParVector *u_edge;
27
28 HYPRE_Int max_iter = maxwell_data-> max_iter;
29 HYPRE_Real tol = maxwell_data-> tol;
30 HYPRE_Int rel_change = maxwell_data-> rel_change;
31 HYPRE_Int zero_guess = maxwell_data-> zero_guess;
32 HYPRE_Int npre_relax = maxwell_data-> num_pre_relax;
33 HYPRE_Int npost_relax = maxwell_data-> num_post_relax;
34
35 hypre_ParCSRMatrix **Ann_l = maxwell_data-> Ann_l;
36 hypre_ParCSRMatrix **Pn_l = maxwell_data-> Pn_l;
37 hypre_ParCSRMatrix **RnT_l = maxwell_data-> RnT_l;
38 hypre_ParVector **bn_l = maxwell_data-> bn_l;
39 hypre_ParVector **xn_l = maxwell_data-> xn_l;
40 hypre_ParVector **resn_l = maxwell_data-> resn_l;
41 hypre_ParVector **en_l = maxwell_data-> en_l;
42 hypre_ParVector **nVtemp_l = maxwell_data-> nVtemp_l;
43 hypre_ParVector **nVtemp2_l = maxwell_data-> nVtemp2_l;
44 HYPRE_Int **nCF_marker_l = maxwell_data-> nCF_marker_l;
45 HYPRE_Real *nrelax_weight= maxwell_data-> nrelax_weight;
46 HYPRE_Real *nomega = maxwell_data-> nomega;
47 HYPRE_Int nrelax_type = maxwell_data-> nrelax_type;
48 HYPRE_Int node_numlevs = maxwell_data-> node_numlevels;
49
50 hypre_ParCSRMatrix *Tgrad = maxwell_data-> Tgrad;
51 hypre_ParCSRMatrix *T_transpose = maxwell_data-> T_transpose;
52
53 hypre_ParCSRMatrix **Aen_l = maxwell_data-> Aen_l;
54 HYPRE_Int en_numlevs = maxwell_data-> en_numlevels;
55
56 hypre_ParCSRMatrix **Aee_l = maxwell_data-> Aee_l;
57 hypre_IJMatrix **Pe_l = maxwell_data-> Pe_l;
58 hypre_IJMatrix **ReT_l = maxwell_data-> ReT_l;
59 hypre_ParVector **be_l = maxwell_data-> be_l;
60 hypre_ParVector **xe_l = maxwell_data-> xe_l;
61 hypre_ParVector **rese_l = maxwell_data-> rese_l;
62 hypre_ParVector **ee_l = maxwell_data-> ee_l;
63 hypre_ParVector **eVtemp_l = maxwell_data-> eVtemp_l;
64 hypre_ParVector **eVtemp2_l = maxwell_data-> eVtemp2_l;
65 HYPRE_Int **eCF_marker_l = maxwell_data-> eCF_marker_l;
66 HYPRE_Real *erelax_weight= maxwell_data-> erelax_weight;
67 HYPRE_Real *eomega = maxwell_data-> eomega;
68 HYPRE_Int erelax_type = maxwell_data-> erelax_type;
69 HYPRE_Int edge_numlevs = maxwell_data-> edge_numlevels;
70
71 HYPRE_Int **BdryRanks_l = maxwell_data-> BdryRanks_l;
72 HYPRE_Int *BdryRanksCnts_l= maxwell_data-> BdryRanksCnts_l;
73
74 HYPRE_Int logging = maxwell_data-> logging;
75 HYPRE_Real *norms = maxwell_data-> norms;
76 HYPRE_Real *rel_norms = maxwell_data-> rel_norms;
77
78 HYPRE_Int relax_local, cycle_param;
79
80 HYPRE_Real b_dot_b = 0, r_dot_r, eps = 0;
81 HYPRE_Real e_dot_e = 0, x_dot_x = 1;
82
83 HYPRE_Int i, j;
84 HYPRE_Int level;
85
86 /* added for the relaxation routines */
87 hypre_ParVector *ze = NULL;
88
89 #if !defined(HYPRE_USING_CUDA) && !defined(HYPRE_USING_HIP)
90 /* GPU impl. needs ze */
91 if (hypre_NumThreads() > 1)
92 #endif
93 {
94 /* Aee is always bigger than Ann */
95
96 ze = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(Aee_l[0]),
97 hypre_ParCSRMatrixGlobalNumRows(Aee_l[0]),
98 hypre_ParCSRMatrixRowStarts(Aee_l[0]));
99 hypre_ParVectorInitialize(ze);
100 }
101
102 hypre_BeginTiming(maxwell_data-> time_index);
103
104 hypre_SStructVectorConvert(f, &f_edge);
105 hypre_SStructVectorConvert(u, &u_edge);
106 hypre_ParVectorZeroBCValues(f_edge, BdryRanks_l[0], BdryRanksCnts_l[0]);
107 hypre_ParVectorZeroBCValues(u_edge, BdryRanks_l[0], BdryRanksCnts_l[0]);
108 be_l[0]= f_edge;
109 xe_l[0]= u_edge;
110
111 /* the nodal fine vectors: bn= T'*be, xn= 0. */
112 hypre_ParCSRMatrixMatvec(1.0, T_transpose, f_edge, 0.0, bn_l[0]);
113 hypre_ParVectorSetConstantValues(xn_l[0], 0.0);
114
115 relax_local= 0;
116 cycle_param= 0;
117
118 (maxwell_data-> num_iterations) = 0;
119 /* if max_iter is zero, return */
120 if (max_iter == 0)
121 {
122 /* if using a zero initial guess, return zero */
123 if (zero_guess)
124 {
125 hypre_ParVectorSetConstantValues(xe_l[0], 0.0);
126 }
127
128 hypre_EndTiming(maxwell_data -> time_index);
129
130 return hypre_error_flag;
131 }
132
133 /* part of convergence check */
134 if (tol > 0.0)
135 {
136 /* eps = (tol^2) */
137 b_dot_b= hypre_ParVectorInnerProd(be_l[0], be_l[0]);
138 eps = tol*tol;
139
140 /* if rhs is zero, return a zero solution */
141 if (b_dot_b == 0.0)
142 {
143 hypre_ParVectorSetConstantValues(xe_l[0], 0.0);
144 if (logging > 0)
145 {
146 norms[0] = 0.0;
147 rel_norms[0] = 0.0;
148 }
149
150 hypre_EndTiming(maxwell_data -> time_index);
151
152 return hypre_error_flag;
153 }
154 }
155
156 /*-----------------------------------------------------
157 * Do V-cycles:
158 * For each index l, "fine" = (l-1), "coarse" = l
159 * down cycle:
160 * a) smooth nodes (Ann)
161 * b) update edge residual (Ane)
162 * c) smooth edges (Aee)
163 * d) restrict updated node and edge residuals
164 * up cycle:
165 * a) interpolate node and edges separately
166 * a) smooth nodes
167 * b) update edge residual
168 * c) smooth edges
169 *
170 * solution update:
171 * edge_sol= edge_sol + T*node_sol
172 *-----------------------------------------------------*/
173 for (i = 0; i < max_iter; i++)
174 {
175 /* fine grid pre_relaxation */
176 for (j = 0; j < npre_relax; j++)
177 {
178 hypre_ParVectorCopy(bn_l[0], nVtemp_l[0]);
179 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[0], xe_l[0],
180 1.0, nVtemp_l[0]);
181
182 hypre_BoomerAMGRelaxIF(Ann_l[0],
183 nVtemp_l[0],
184 nCF_marker_l[0],
185 nrelax_type,
186 relax_local,
187 cycle_param,
188 nrelax_weight[0],
189 nomega[0],
190 NULL,
191 xn_l[0],
192 nVtemp2_l[0],
193 ze);
194
195 /* update edge right-hand fe_l= fe_l-Aen_l*xn_l[0] */
196 hypre_ParVectorCopy(be_l[0], eVtemp_l[0]);
197 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[0], xn_l[0],
198 1.0, eVtemp_l[0]);
199 hypre_ParVectorZeroBCValues(eVtemp_l[0], BdryRanks_l[0],
200 BdryRanksCnts_l[0]);
201
202 hypre_BoomerAMGRelaxIF(Aee_l[0],
203 eVtemp_l[0],
204 eCF_marker_l[0],
205 erelax_type,
206 relax_local,
207 cycle_param,
208 erelax_weight[0],
209 eomega[0],
210 NULL,
211 xe_l[0],
212 eVtemp2_l[0],
213 ze);
214 } /* for (j = 0; j < npre_relax; j++) */
215
216 /* compute fine grid residual. Note the edge residual of
217 the block system is the residual of the actual edge equations
218 itself. */
219 hypre_ParVectorCopy(bn_l[0], resn_l[0]);
220 hypre_ParCSRMatrixMatvec(-1.0, Ann_l[0], xn_l[0], 1.0, resn_l[0]);
221 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[0], xe_l[0], 1.0, resn_l[0]);
222
223 hypre_ParVectorCopy(be_l[0], rese_l[0]);
224 hypre_ParCSRMatrixMatvec(-1.0, Aee_l[0], xe_l[0], 1.0, rese_l[0]);
225 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[0], xn_l[0], 1.0, rese_l[0]);
226 hypre_ParVectorZeroBCValues(rese_l[0], BdryRanks_l[0], BdryRanksCnts_l[0]);
227
228 /* convergence check */
229 if (tol > 0.0)
230 {
231 r_dot_r= hypre_ParVectorInnerProd(rese_l[0], rese_l[0]);
232
233 if (logging > 0)
234 {
235 norms[i] = sqrt(r_dot_r);
236 if (b_dot_b > 0)
237 {
238 rel_norms[i] = sqrt(r_dot_r/b_dot_b);
239 }
240 else
241 {
242 rel_norms[i] = 0.0;
243 }
244 }
245
246 /* always do at least 1 V-cycle */
247 if ((r_dot_r/b_dot_b < eps) && (i > 0))
248 {
249 if (rel_change)
250 {
251 if ((e_dot_e/x_dot_x) < eps)
252 {
253 break;
254 }
255 }
256 else
257 {
258 break;
259 }
260 }
261 }
262
263 if (en_numlevs > 1)
264 {
265 hypre_ParCSRMatrixMatvecT(1.0, RnT_l[0], resn_l[0], 0.0,
266 bn_l[1]);
267
268 hypre_ParCSRMatrixMatvecT(1.0,
269 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(ReT_l[0]),
270 rese_l[0], 0.0, be_l[1]);
271
272 hypre_ParVectorZeroBCValues(be_l[1], BdryRanks_l[1],
273 BdryRanksCnts_l[1]);
274
275 /* zero off initial guess for the next level */
276 hypre_ParVectorSetConstantValues(xn_l[1], 0.0);
277 hypre_ParVectorSetConstantValues(xe_l[1], 0.0);
278
279 } /* if (en_numlevs > 1) */
280
281 for (level = 1; level<= en_numlevs-2; level++)
282 {
283 /*-----------------------------------------------
284 * Down cycle
285 *-----------------------------------------------*/
286 for (j = 0; j < npre_relax; j++)
287 {
288 hypre_ParVectorCopy(bn_l[level], nVtemp_l[level]);
289 if (j)
290 {
291 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[level],
292 xe_l[level], 1.0, nVtemp_l[level]);
293 }
294 hypre_BoomerAMGRelaxIF(Ann_l[level],
295 nVtemp_l[level],
296 nCF_marker_l[level],
297 nrelax_type,
298 relax_local,
299 cycle_param,
300 nrelax_weight[level],
301 nomega[level],
302 NULL,
303 xn_l[level],
304 nVtemp2_l[level],
305 ze);
306
307 /* update edge right-hand fe_l= fe_l-Aen_l*xn_l[level] */
308 hypre_ParVectorCopy(be_l[level], eVtemp_l[level]);
309 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[level],
310 xn_l[level], 1.0, eVtemp_l[level]);
311 hypre_ParVectorZeroBCValues(eVtemp_l[level], BdryRanks_l[level],
312 BdryRanksCnts_l[level]);
313
314 hypre_BoomerAMGRelaxIF(Aee_l[level],
315 eVtemp_l[level],
316 eCF_marker_l[level],
317 erelax_type,
318 relax_local,
319 cycle_param,
320 erelax_weight[level],
321 eomega[level],
322 NULL,
323 xe_l[level],
324 eVtemp2_l[level],
325 ze);
326 } /*for (j = 0; j < npre_relax; j++) */
327
328 /* compute residuals */
329 hypre_ParVectorCopy(bn_l[level], resn_l[level]);
330 hypre_ParCSRMatrixMatvec(-1.0, Ann_l[level], xn_l[level],
331 1.0, resn_l[level]);
332
333 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[level], xe_l[level],
334 1.0, resn_l[level]);
335
336 hypre_ParVectorCopy(be_l[level], rese_l[level]);
337 hypre_ParCSRMatrixMatvec(-1.0, Aee_l[level], xe_l[level],
338 1.0, rese_l[level]);
339 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[level], xn_l[level],
340 1.0, rese_l[level]);
341 hypre_ParVectorZeroBCValues(rese_l[level], BdryRanks_l[level],
342 BdryRanksCnts_l[level]);
343
344 /* restrict residuals */
345 hypre_ParCSRMatrixMatvecT(1.0, RnT_l[level], resn_l[level],
346 0.0, bn_l[level+1]);
347
348 hypre_ParCSRMatrixMatvecT(1.0,
349 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(ReT_l[level]),
350 rese_l[level], 0.0, be_l[level+1]);
351
352 hypre_ParVectorZeroBCValues(be_l[level+1], BdryRanks_l[level+1],
353 BdryRanksCnts_l[level+1]);
354
355 /* zero off initial guess for the next level */
356 hypre_ParVectorSetConstantValues(xn_l[level+1], 0.0);
357 hypre_ParVectorSetConstantValues(xe_l[level+1], 0.0);
358
359 } /* for (level = 0; level<= en_numlevels-2; level++) */
360
361 /*----------------------------------------------------------------
362 * For the lowest edge-node level, solve using relaxation or
363 * cycling down if there are more than en_numlevels levels for
364 * one of the node or edge dofs.
365 *----------------------------------------------------------------*/
366 level = en_numlevs-1;
367
368 /* npre_relax if not the coarsest level. Otherwise, relax once.*/
369 if ( (en_numlevs != edge_numlevs)
370 || (en_numlevs != node_numlevs) )
371 {
372 for (j = 0; j < npre_relax; j++)
373 {
374 hypre_ParVectorCopy(bn_l[level], nVtemp_l[level]);
375 if (j)
376 {
377 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[level],
378 xe_l[level], 1.0, nVtemp_l[level]);
379 }
380 hypre_BoomerAMGRelaxIF(Ann_l[level],
381 nVtemp_l[level],
382 nCF_marker_l[level],
383 nrelax_type,
384 relax_local,
385 cycle_param,
386 nrelax_weight[level],
387 nomega[level],
388 NULL,
389 xn_l[level],
390 nVtemp2_l[level],
391 ze);
392
393 /* update edge right-hand fe_l= fe_l-Aen_l*xn_l[level] */
394 hypre_ParVectorCopy(be_l[level], eVtemp_l[level]);
395 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[level],
396 xn_l[level], 1.0, eVtemp_l[level]);
397
398 hypre_ParVectorZeroBCValues(eVtemp_l[level], BdryRanks_l[level],
399 BdryRanksCnts_l[level]);
400
401 hypre_BoomerAMGRelaxIF(Aee_l[level],
402 eVtemp_l[level],
403 eCF_marker_l[level],
404 erelax_type,
405 relax_local,
406 cycle_param,
407 erelax_weight[level],
408 eomega[level],
409 NULL,
410 xe_l[level],
411 eVtemp2_l[level],
412 ze);
413 } /*for (j = 0; j < npre_relax; j++) */
414 } /* if ( (en_numlevs != edge_numlevs) */
415
416 else
417 {
418 hypre_BoomerAMGRelaxIF(Ann_l[level],
419 bn_l[level],
420 nCF_marker_l[level],
421 nrelax_type,
422 relax_local,
423 cycle_param,
424 nrelax_weight[level],
425 nomega[level],
426 NULL,
427 xn_l[level],
428 nVtemp2_l[level],
429 ze);
430
431 hypre_ParVectorCopy(be_l[level], eVtemp_l[level]);
432 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[level], xn_l[level],
433 1.0, eVtemp_l[level]);
434
435 hypre_ParVectorZeroBCValues(eVtemp_l[level], BdryRanks_l[level],
436 BdryRanksCnts_l[level]);
437
438 hypre_BoomerAMGRelaxIF(Aee_l[level],
439 eVtemp_l[level],
440 eCF_marker_l[level],
441 erelax_type,
442 relax_local,
443 cycle_param,
444 erelax_weight[level],
445 eomega[level],
446 NULL,
447 xe_l[level],
448 eVtemp2_l[level],
449 ze);
450 }
451
452 /* Continue down the edge hierarchy if more edge levels. */
453 if (edge_numlevs > en_numlevs)
454 {
455 hypre_ParVectorCopy(be_l[level], rese_l[level]);
456 hypre_ParCSRMatrixMatvec(-1.0, Aee_l[level], xe_l[level], 1.0,
457 rese_l[level]);
458 hypre_ParCSRMatrixMatvecT(1.0,
459 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(ReT_l[level]),
460 rese_l[level], 0.0, be_l[level+1]);
461 hypre_ParVectorZeroBCValues(be_l[level+1], BdryRanks_l[level+1],
462 BdryRanksCnts_l[level+1]);
463
464 hypre_ParVectorSetConstantValues(xe_l[level+1], 0.0);
465
466 for (level = en_numlevs; level<= edge_numlevs-2; level++)
467 {
468 for (j = 0; j < npre_relax; j++)
469 {
470 hypre_BoomerAMGRelaxIF(Aee_l[level],
471 be_l[level],
472 eCF_marker_l[level],
473 erelax_type,
474 relax_local,
475 cycle_param,
476 erelax_weight[level],
477 eomega[level],
478 NULL,
479 xe_l[level],
480 eVtemp2_l[level],
481 ze);
482 }
483
484 /* compute residuals and restrict */
485 hypre_ParVectorCopy(be_l[level], rese_l[level]);
486 hypre_ParCSRMatrixMatvec(-1.0, Aee_l[level], xe_l[level],
487 1.0, rese_l[level]);
488 hypre_ParCSRMatrixMatvecT(1.0,
489 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(ReT_l[level]),
490 rese_l[level], 0.0, be_l[level+1]);
491 hypre_ParVectorZeroBCValues(be_l[level+1], BdryRanks_l[level+1],
492 BdryRanksCnts_l[level+1]);
493
494 hypre_ParVectorSetConstantValues(xe_l[level+1], 0.0);
495 } /* for (level = en_numlevs; level< edge_numlevs-2; level++) */
496
497 /* coarsest relaxation */
498 level = edge_numlevs-1;
499 hypre_BoomerAMGRelaxIF(Aee_l[level],
500 be_l[level],
501 eCF_marker_l[level],
502 erelax_type,
503 relax_local,
504 cycle_param,
505 erelax_weight[level],
506 eomega[level],
507 NULL,
508 xe_l[level],
509 eVtemp2_l[level],
510 ze);
511 } /* if (edge_numlevs > en_numlevs) */
512
513 /*-----------------------------------------------------------
514 * node hierarchy has more levels than the edge hierarchy:
515 * continue to march down the node hierarchy
516 *-----------------------------------------------------------*/
517 else if (node_numlevs > en_numlevs)
518 {
519 hypre_ParVectorCopy(bn_l[level], resn_l[level]);
520 hypre_ParCSRMatrixMatvec(-1.0, Ann_l[level], xn_l[level], 1.0,
521 resn_l[level]);
522 hypre_ParCSRMatrixMatvecT(1.0,
523 (hypre_ParCSRMatrix *) RnT_l[level],
524 resn_l[level], 0.0, bn_l[level+1]);
525
526 hypre_ParVectorSetConstantValues(xn_l[level+1], 0.0);
527
528 for (level = en_numlevs; level<= node_numlevs-2; level++)
529 {
530 for (j = 0; j < npre_relax; j++)
531 {
532 hypre_BoomerAMGRelaxIF(Ann_l[level],
533 bn_l[level],
534 nCF_marker_l[level],
535 nrelax_type,
536 relax_local,
537 cycle_param,
538 nrelax_weight[level],
539 nomega[level],
540 NULL,
541 xn_l[level],
542 nVtemp2_l[level],
543 ze);
544 }
545
546 /* compute residuals and restrict */
547 hypre_ParVectorCopy(bn_l[level], resn_l[level]);
548 hypre_ParCSRMatrixMatvec(-1.0, Ann_l[level], xn_l[level],
549 1.0, resn_l[level]);
550 hypre_ParCSRMatrixMatvecT(1.0, RnT_l[level], resn_l[level],
551 0.0, bn_l[level+1]);
552
553 hypre_ParVectorSetConstantValues(xn_l[level+1], 0.0);
554 } /* for (level = en_numlevs; level<= node_numlevs-2; level++) */
555
556 /* coarsest relaxation */
557 level = node_numlevs-1;
558 hypre_BoomerAMGRelaxIF(Ann_l[level],
559 bn_l[level],
560 nCF_marker_l[level],
561 nrelax_type,
562 relax_local,
563 cycle_param,
564 nrelax_weight[level],
565 nomega[level],
566 NULL,
567 xn_l[level],
568 nVtemp2_l[level],
569 ze);
570 } /* else if (node_numlevs > en_numlevs) */
571
572 /*---------------------------------------------------------------------
573 * Up cycle. First the extra hierarchy levels. Notice we relax on
574 * the coarsest en_numlevel.
575 *---------------------------------------------------------------------*/
576 if (edge_numlevs > en_numlevs)
577 {
578 for (level = (edge_numlevs - 2); level>= en_numlevs-1; level--)
579 {
580 hypre_ParCSRMatrixMatvec(1.0,
581 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(Pe_l[level]),
582 xe_l[level+1], 0.0, ee_l[level]);
583 hypre_ParVectorZeroBCValues(ee_l[level], BdryRanks_l[level],
584 BdryRanksCnts_l[level]);
585 hypre_ParVectorAxpy(1.0, ee_l[level], xe_l[level]);
586
587 /* post smooth */
588 for (j = 0; j < npost_relax; j++)
589 {
590 hypre_BoomerAMGRelaxIF(Aee_l[level],
591 be_l[level],
592 eCF_marker_l[level],
593 erelax_type,
594 relax_local,
595 cycle_param,
596 erelax_weight[level],
597 eomega[level],
598 NULL,
599 xe_l[level],
600 eVtemp2_l[level],
601 ze);
602 }
603
604 } /* for (level = (edge_numlevs - 2); level>= en_numlevs; level--) */
605 } /* if (edge_numlevs > en_numlevs) */
606
607 else if (node_numlevs > en_numlevs)
608 {
609 for (level = (node_numlevs - 2); level>= en_numlevs-1; level--)
610 {
611 hypre_ParCSRMatrixMatvec(1.0, Pn_l[level], xn_l[level+1], 0.0,
612 en_l[level]);
613 hypre_ParVectorAxpy(1.0, en_l[level], xn_l[level]);
614
615 /* post smooth */
616 for (j = 0; j < npost_relax; j++)
617 {
618 hypre_BoomerAMGRelaxIF(Ann_l[level],
619 bn_l[level],
620 nCF_marker_l[level],
621 nrelax_type,
622 relax_local,
623 cycle_param,
624 nrelax_weight[level],
625 nomega[level],
626 NULL,
627 xn_l[level],
628 nVtemp2_l[level],
629 ze);
630 }
631
632 } /* for (level = (node_numlevs - 2); level>= en_numlevs; level--) */
633 } /* else if (node_numlevs > en_numlevs) */
634
635 /*---------------------------------------------------------------------
636 * Cycle up the common levels.
637 *---------------------------------------------------------------------*/
638 for (level = (en_numlevs - 2); level>= 1; level--)
639 {
640 hypre_ParCSRMatrixMatvec(1.0, Pn_l[level], xn_l[level+1], 0.0,
641 en_l[level]);
642 hypre_ParVectorAxpy(1.0, en_l[level], xn_l[level]);
643
644 hypre_ParCSRMatrixMatvec(1.0,
645 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(Pe_l[level]),
646 xe_l[level+1], 0.0, ee_l[level]);
647 hypre_ParVectorZeroBCValues(ee_l[level], BdryRanks_l[level],
648 BdryRanksCnts_l[level]);
649 hypre_ParVectorAxpy(1.0, ee_l[level], xe_l[level]);
650
651 /* post smooth */
652 for (j = 0; j < npost_relax; j++)
653 {
654 hypre_ParVectorCopy(bn_l[level], nVtemp_l[level]);
655 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[level], xe_l[level],
656 1.0, nVtemp_l[level]);
657 hypre_BoomerAMGRelaxIF(Ann_l[level],
658 nVtemp_l[level],
659 nCF_marker_l[level],
660 nrelax_type,
661 relax_local,
662 cycle_param,
663 nrelax_weight[level],
664 nomega[level],
665 NULL,
666 xn_l[level],
667 nVtemp_l[level],
668 ze);
669
670 hypre_ParVectorCopy(be_l[level], eVtemp_l[level]);
671 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[level], xn_l[level],
672 1.0, eVtemp_l[level]);
673 hypre_ParVectorZeroBCValues(eVtemp_l[level], BdryRanks_l[level],
674 BdryRanksCnts_l[level]);
675
676 hypre_BoomerAMGRelaxIF(Aee_l[level],
677 eVtemp_l[level],
678 eCF_marker_l[level],
679 erelax_type,
680 relax_local,
681 cycle_param,
682 erelax_weight[level],
683 eomega[level],
684 NULL,
685 xe_l[level],
686 eVtemp2_l[level],
687 ze);
688 }
689
690 } /* for (level = (en_numlevs - 2); level>= 1; level--) */
691
692 /* interpolate error and correct on finest grids */
693 hypre_ParCSRMatrixMatvec(1.0, Pn_l[0], xn_l[1], 0.0, en_l[0]);
694 hypre_ParVectorAxpy(1.0, en_l[0], xn_l[0]);
695
696 hypre_ParCSRMatrixMatvec(1.0,
697 (hypre_ParCSRMatrix *) hypre_IJMatrixObject(Pe_l[0]),
698 xe_l[1], 0.0, ee_l[0]);
699 hypre_ParVectorZeroBCValues(ee_l[0], BdryRanks_l[0],
700 BdryRanksCnts_l[0]);
701 hypre_ParVectorAxpy(1.0, ee_l[0], xe_l[0]);
702
703 /* part of convergence check. Will assume that if en_numlevels= 1,
704 then so would edge_numlevels and node_numlevels. Otherwise,
705 we measure the error of xe_l[0] + T*xn_l[0]. */
706 if ((tol > 0.0) && (rel_change))
707 {
708 if (en_numlevs > 1)
709 {
710 hypre_ParCSRMatrixMatvec(1.0, Tgrad, en_l[0], 1.0,
711 ee_l[0]);
712 hypre_ParVectorZeroBCValues(ee_l[0], BdryRanks_l[0],
713 BdryRanksCnts_l[0]);
714 e_dot_e= hypre_ParVectorInnerProd(ee_l[0], ee_l[0]);
715
716 hypre_ParVectorCopy(xe_l[0], eVtemp_l[0]);
717 hypre_ParCSRMatrixMatvec(1.0, Tgrad, xn_l[0], 1.0,
718 eVtemp_l[0]);
719 hypre_ParVectorZeroBCValues(eVtemp_l[0], BdryRanks_l[0],
720 BdryRanksCnts_l[0]);
721 x_dot_x= hypre_ParVectorInnerProd(eVtemp_l[0], eVtemp_l[0]);
722 }
723 else
724 {
725 e_dot_e = 0.0;
726 x_dot_x = 1.0;
727 }
728 }
729
730 /* check nodal convergence */
731
732 for (j = 0; j < npost_relax; j++)
733 {
734 hypre_ParVectorCopy(bn_l[0], nVtemp_l[0]);
735 hypre_ParCSRMatrixMatvecT(-1.0, Aen_l[0], xe_l[0],
736 1.0, nVtemp_l[0]);
737 hypre_BoomerAMGRelaxIF(Ann_l[0],
738 nVtemp_l[0],
739 nCF_marker_l[0],
740 nrelax_type,
741 relax_local,
742 cycle_param,
743 nrelax_weight[0],
744 nomega[0],
745 NULL,
746 xn_l[0],
747 nVtemp2_l[0],
748 ze);
749
750 hypre_ParVectorCopy(be_l[0], eVtemp_l[0]);
751 hypre_ParCSRMatrixMatvec(-1.0, Aen_l[0], xn_l[0], 1.0,
752 eVtemp_l[0]);
753 hypre_ParVectorZeroBCValues(eVtemp_l[0], BdryRanks_l[0],
754 BdryRanksCnts_l[0]);
755
756 hypre_BoomerAMGRelaxIF(Aee_l[0],
757 eVtemp_l[0],
758 eCF_marker_l[0],
759 erelax_type,
760 relax_local,
761 cycle_param,
762 erelax_weight[0],
763 eomega[0],
764 NULL,
765 xe_l[0],
766 eVtemp2_l[0],
767 ze);
768 } /* for (j = 0; j < npost_relax; j++) */
769
770 (maxwell_data -> num_iterations) = (i + 1);
771 }
772
773 /* add the gradient solution component to u_edge */
774 hypre_ParCSRMatrixMatvec(1.0, Tgrad, xn_l[0], 1.0, u_edge);
775 hypre_ParVectorZeroBCValues(u_edge, BdryRanks_l[0], BdryRanksCnts_l[0]);
776
777 hypre_EndTiming(maxwell_data -> time_index);
778 hypre_ParVectorDestroy(ze);
779
780 return hypre_error_flag;
781 }
782