source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/gen_redcs_mat.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: 15.9 KB
Line 
1#include "parcsr_ls.h"
2#include "parcsr_mv.h"
3#include "seq_mv.h"
4#include "par_amg.h"
5
6
7/* here we have the sequential setup and solve - called from the
8 * parallel one - for the coarser levels */
9
10int hypre_seqAMGSetup( hypre_ParAMGData *amg_data,
11 int p_level,
12 int coarse_threshold)
13
14
15{
16
17 /* Par Data Structure variables */
18 hypre_ParCSRMatrix **Par_A_array = hypre_ParAMGDataAArray(amg_data);
19
20 MPI_Comm comm = hypre_ParCSRMatrixComm(Par_A_array[0]);
21 MPI_Comm new_comm, seq_comm;
22
23 hypre_ParCSRMatrix *A_seq = NULL;
24 hypre_CSRMatrix *A_seq_diag;
25 hypre_CSRMatrix *A_seq_offd;
26 hypre_ParVector *F_seq = NULL;
27 hypre_ParVector *U_seq = NULL;
28
29 hypre_ParCSRMatrix *A;
30
31 int **dof_func_array;
32 int num_procs, my_id;
33
34 int not_finished_coarsening;
35 int level;
36 int redundant;
37
38 HYPRE_Solver coarse_solver;
39
40 /* misc */
41 dof_func_array = hypre_ParAMGDataDofFuncArray(amg_data);
42 redundant = hypre_ParAMGDataRedundant(amg_data);
43
44 /*MPI Stuff */
45 MPI_Comm_size(comm, &num_procs);
46
47 /*initial */
48 level = p_level;
49
50 not_finished_coarsening = 1;
51
52 /* convert A at this level to sequential */
53 A = Par_A_array[level];
54
55 {
56 double *A_seq_data = NULL;
57 int *A_seq_i = NULL;
58 int *A_seq_offd_i = NULL;
59 int *A_seq_j = NULL;
60
61 double *A_tmp_data = NULL;
62 int *A_tmp_i = NULL;
63 HYPRE_BigInt *A_tmp_j = NULL;
64
65 int *info = NULL;
66 int *displs = NULL;
67 int *displs2 = NULL;
68 int i, j, size, num_nonzeros, total_nnz, cnt;
69
70 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
71 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
72 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
73 int *A_diag_i = hypre_CSRMatrixI(A_diag);
74 int *A_offd_i = hypre_CSRMatrixI(A_offd);
75 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
76 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
77 double *A_diag_data = hypre_CSRMatrixData(A_diag);
78 double *A_offd_data = hypre_CSRMatrixData(A_offd);
79 int num_rows = hypre_CSRMatrixNumRows(A_diag);
80 HYPRE_BigInt first_row_index = hypre_ParCSRMatrixFirstRowIndex(A);
81 int new_num_procs;
82 HYPRE_BigInt *row_starts;
83
84 hypre_GenerateSubComm(comm, num_rows, &new_comm);
85
86 if (num_rows)
87 {
88 hypre_ParAMGDataParticipate(amg_data) = 1;
89 MPI_Comm_size(new_comm, &new_num_procs);
90 MPI_Comm_rank(new_comm, &my_id);
91 info = hypre_CTAlloc(int, new_num_procs);
92
93 if (redundant)
94 MPI_Allgather(&num_rows, 1, MPI_INT, info, 1, MPI_INT, new_comm);
95 else
96 MPI_Gather(&num_rows, 1, MPI_INT, info, 1, MPI_INT, 0, new_comm);
97
98 /* alloc space in seq data structure only for participating procs*/
99 if (redundant || my_id == 0)
100 {
101 HYPRE_BoomerAMGCreate(&coarse_solver);
102 HYPRE_BoomerAMGSetMaxRowSum(coarse_solver,
103 hypre_ParAMGDataMaxRowSum(amg_data));
104 HYPRE_BoomerAMGSetStrongThreshold(coarse_solver,
105 hypre_ParAMGDataStrongThreshold(amg_data));
106 HYPRE_BoomerAMGSetCoarsenType(coarse_solver,
107 hypre_ParAMGDataCoarsenType(amg_data));
108 HYPRE_BoomerAMGSetInterpType(coarse_solver,
109 hypre_ParAMGDataInterpType(amg_data));
110 HYPRE_BoomerAMGSetTruncFactor(coarse_solver,
111 hypre_ParAMGDataTruncFactor(amg_data));
112 HYPRE_BoomerAMGSetPMaxElmts(coarse_solver,
113 hypre_ParAMGDataPMaxElmts(amg_data));
114 HYPRE_BoomerAMGSetGridRelaxType(coarse_solver,
115 hypre_ParAMGDataGridRelaxType(amg_data));
116 HYPRE_BoomerAMGSetGridRelaxPoints(coarse_solver,
117 hypre_ParAMGDataGridRelaxPoints(amg_data));
118 HYPRE_BoomerAMGSetRelaxOrder(coarse_solver,
119 hypre_ParAMGDataRelaxOrder(amg_data));
120 HYPRE_BoomerAMGSetRelaxWeight(coarse_solver,
121 &hypre_ParAMGDataRelaxWeight(amg_data)[p_level]);
122 HYPRE_BoomerAMGSetNumGridSweeps(coarse_solver,
123 hypre_ParAMGDataNumGridSweeps(amg_data));
124 HYPRE_BoomerAMGSetNumFunctions(coarse_solver,
125 hypre_ParAMGDataNumFunctions(amg_data));
126 HYPRE_BoomerAMGSetMaxIter(coarse_solver, 1);
127 HYPRE_BoomerAMGSetTol(coarse_solver, 0);
128 }
129
130 /* Create CSR Matrix, will be Diag part of new matrix */
131 A_tmp_i = hypre_CTAlloc(int, num_rows+1);
132
133 A_tmp_i[0] = 0;
134 for (i=1; i < num_rows+1; i++)
135 A_tmp_i[i] = A_diag_i[i]-A_diag_i[i-1]+A_offd_i[i]-A_offd_i[i-1];
136
137 num_nonzeros = A_offd_i[num_rows]+A_diag_i[num_rows];
138
139 A_tmp_j = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros);
140 A_tmp_data = hypre_CTAlloc(double, num_nonzeros);
141
142 cnt = 0;
143 for (i=0; i < num_rows; i++)
144 {
145 for (j=A_diag_i[i]; j < A_diag_i[i+1]; j++)
146 {
147 A_tmp_j[cnt] = (HYPRE_BigInt) A_diag_j[j]+first_row_index;
148 A_tmp_data[cnt++] = A_diag_data[j];
149 }
150 for (j=A_offd_i[i]; j < A_offd_i[i+1]; j++)
151 {
152 A_tmp_j[cnt] = col_map_offd[A_offd_j[j]];
153 A_tmp_data[cnt++] = A_offd_data[j];
154 }
155 }
156
157 displs = hypre_CTAlloc(int, new_num_procs+1);
158 displs[0] = 0;
159 for (i=1; i < new_num_procs+1; i++)
160 displs[i] = displs[i-1]+info[i-1];
161 size = displs[new_num_procs];
162
163 if (redundant || my_id == 0)
164 {
165 A_seq_i = hypre_CTAlloc(int, size+1);
166 A_seq_offd_i = hypre_CTAlloc(int, size+1);
167 }
168
169 if (redundant)
170 MPI_Allgatherv ( &A_tmp_i[1], num_rows, MPI_INT, &A_seq_i[1], info,
171 displs, MPI_INT, new_comm );
172 else
173 MPI_Gatherv ( &A_tmp_i[1], num_rows, MPI_INT, &A_seq_i[1], info,
174 displs, MPI_INT, 0, new_comm );
175
176 if (redundant || my_id == 0)
177 {
178 displs2 = hypre_CTAlloc(int, new_num_procs+1);
179
180 A_seq_i[0] = 0;
181 displs2[0] = 0;
182 for (j=1; j < displs[1]; j++)
183 A_seq_i[j] = A_seq_i[j]+A_seq_i[j-1];
184 for (i=1; i < new_num_procs; i++)
185 {
186 for (j=displs[i]; j < displs[i+1]; j++)
187 {
188 A_seq_i[j] = A_seq_i[j]+A_seq_i[j-1];
189 }
190 }
191 A_seq_i[size] = A_seq_i[size]+A_seq_i[size-1];
192 displs2[new_num_procs] = A_seq_i[size];
193 for (i=1; i < new_num_procs+1; i++)
194 {
195 displs2[i] = A_seq_i[displs[i]];
196 info[i-1] = displs2[i] - displs2[i-1];
197 }
198
199 total_nnz = displs2[new_num_procs];
200 A_seq_j = hypre_CTAlloc(int, total_nnz);
201 A_seq_data = hypre_CTAlloc(double, total_nnz);
202 }
203 if (redundant)
204 {
205 MPI_Allgatherv ( A_tmp_j, num_nonzeros, MPI_HYPRE_BIG_INT,
206 A_seq_j, info, displs2,
207 MPI_INT, new_comm );
208
209 MPI_Allgatherv ( A_tmp_data, num_nonzeros, MPI_DOUBLE,
210 A_seq_data, info, displs2,
211 MPI_DOUBLE, new_comm );
212 }
213 else
214 {
215 MPI_Gatherv ( A_tmp_j, num_nonzeros, MPI_HYPRE_BIG_INT,
216 A_seq_j, info, displs2,
217 MPI_INT, 0, new_comm );
218
219 MPI_Gatherv ( A_tmp_data, num_nonzeros, MPI_DOUBLE,
220 A_seq_data, info, displs2,
221 MPI_DOUBLE, 0, new_comm );
222 }
223
224 hypre_TFree(info);
225 hypre_TFree(displs);
226 hypre_TFree(A_tmp_i);
227 hypre_TFree(A_tmp_j);
228 hypre_TFree(A_tmp_data);
229
230 if (redundant || my_id == 0)
231 {
232 hypre_TFree(displs2);
233
234 row_starts = hypre_CTAlloc(HYPRE_BigInt,2);
235 row_starts[0] = 0;
236 row_starts[1] = size;
237
238 /* Create 1 proc communicator */
239 seq_comm = MPI_COMM_SELF;
240
241 A_seq = hypre_ParCSRMatrixCreate(seq_comm,size,size,
242 row_starts, row_starts,
243 0,total_nnz,0);
244
245 A_seq_diag = hypre_ParCSRMatrixDiag(A_seq);
246 A_seq_offd = hypre_ParCSRMatrixOffd(A_seq);
247
248 hypre_CSRMatrixData(A_seq_diag) = A_seq_data;
249 hypre_CSRMatrixI(A_seq_diag) = A_seq_i;
250 hypre_CSRMatrixJ(A_seq_diag) = A_seq_j;
251 hypre_CSRMatrixI(A_seq_offd) = A_seq_offd_i;
252
253 F_seq = hypre_ParVectorCreate(seq_comm, size, row_starts);
254 U_seq = hypre_ParVectorCreate(seq_comm, size, row_starts);
255 hypre_ParVectorOwnsPartitioning(F_seq) = 0;
256 hypre_ParVectorOwnsPartitioning(U_seq) = 0;
257 hypre_ParVectorInitialize(F_seq);
258 hypre_ParVectorInitialize(U_seq);
259
260 hypre_BoomerAMGSetup(coarse_solver,A_seq,F_seq,U_seq);
261
262 hypre_ParAMGDataCoarseSolver(amg_data) = coarse_solver;
263 hypre_ParAMGDataACoarse(amg_data) = A_seq;
264 hypre_ParAMGDataFCoarse(amg_data) = F_seq;
265 hypre_ParAMGDataUCoarse(amg_data) = U_seq;
266 }
267 hypre_ParAMGDataNewComm(amg_data) = new_comm;
268 }
269 }
270 return 0;
271}
272
273
274
275/*--------------------------------------------------------------------------
276 * hypre_seqAMGCycle
277 *--------------------------------------------------------------------------*/
278
279int
280hypre_seqAMGCycle( hypre_ParAMGData *amg_data,
281 int p_level,
282 hypre_ParVector **Par_F_array,
283 hypre_ParVector **Par_U_array )
284{
285
286 hypre_ParVector *Aux_U;
287 hypre_ParVector *Aux_F;
288
289 /* Local variables */
290
291 int Solve_err_flag = 0;
292
293 int n;
294 int i;
295
296 hypre_Vector *u_local;
297 double *u_data;
298
299 HYPRE_BigInt first_index;
300
301 /* Acquire seq data */
302 MPI_Comm new_comm = hypre_ParAMGDataNewComm(amg_data);
303 HYPRE_Solver coarse_solver = hypre_ParAMGDataCoarseSolver(amg_data);
304 hypre_ParCSRMatrix *A_coarse = hypre_ParAMGDataACoarse(amg_data);
305 hypre_ParVector *F_coarse = hypre_ParAMGDataFCoarse(amg_data);
306 hypre_ParVector *U_coarse = hypre_ParAMGDataUCoarse(amg_data);
307 int redundant = hypre_ParAMGDataRedundant(amg_data);
308
309 Aux_U = Par_U_array[p_level];
310 Aux_F = Par_F_array[p_level];
311
312 first_index = hypre_ParVectorFirstIndex(Aux_U);
313 u_local = hypre_ParVectorLocalVector(Aux_U);
314 u_data = hypre_VectorData(u_local);
315 n = hypre_VectorSize(u_local);
316
317
318 /*if (A_coarse)*/
319 if (hypre_ParAMGDataParticipate(amg_data))
320 {
321 double *f_data;
322 hypre_Vector *f_local;
323 hypre_Vector *tmp_vec;
324
325 int nf;
326 int local_info;
327 double *recv_buf;
328 int *displs = NULL;
329 int *info = NULL;
330 int size;
331 int new_num_procs, my_id;
332
333 MPI_Comm_size(new_comm, &new_num_procs);
334 MPI_Comm_rank(new_comm, &my_id);
335
336 f_local = hypre_ParVectorLocalVector(Aux_F);
337 f_data = hypre_VectorData(f_local);
338 nf = hypre_VectorSize(f_local);
339
340 /* first f */
341 info = hypre_CTAlloc(int, new_num_procs);
342 local_info = nf;
343 if (redundant)
344 MPI_Allgather(&local_info, 1, MPI_INT, info, 1, MPI_INT, new_comm);
345 else
346 MPI_Gather(&local_info, 1, MPI_INT, info, 1, MPI_INT, 0, new_comm);
347
348 if (redundant || my_id ==0)
349 {
350 displs = hypre_CTAlloc(int, new_num_procs+1);
351 displs[0] = 0;
352 for (i=1; i < new_num_procs+1; i++)
353 displs[i] = displs[i-1]+info[i-1];
354 size = displs[new_num_procs];
355
356 tmp_vec = hypre_ParVectorLocalVector(F_coarse);
357 recv_buf = hypre_VectorData(tmp_vec);
358 }
359
360 if (redundant)
361 MPI_Allgatherv ( f_data, nf, MPI_DOUBLE,
362 recv_buf, info, displs,
363 MPI_DOUBLE, new_comm );
364 else
365 MPI_Gatherv ( f_data, nf, MPI_DOUBLE,
366 recv_buf, info, displs,
367 MPI_DOUBLE, 0, new_comm );
368
369 if (redundant || my_id ==0)
370 {
371 tmp_vec = hypre_ParVectorLocalVector(U_coarse);
372 recv_buf = hypre_VectorData(tmp_vec);
373 }
374
375 /*then u */
376 if (redundant)
377 {
378 MPI_Allgatherv ( u_data, n, MPI_DOUBLE,
379 recv_buf, info, displs,
380 MPI_DOUBLE, new_comm );
381 hypre_TFree(displs);
382 hypre_TFree(info);
383 }
384 else
385 MPI_Gatherv ( u_data, n, MPI_DOUBLE,
386 recv_buf, info, displs,
387 MPI_DOUBLE, 0, new_comm );
388
389 /* clean up */
390 if (redundant || my_id ==0)
391 {
392 hypre_BoomerAMGSolve(coarse_solver, A_coarse, F_coarse, U_coarse);
393 }
394
395 /*copy my part of U to parallel vector */
396 if (redundant)
397 {
398 double *local_data;
399
400 local_data = hypre_VectorData(hypre_ParVectorLocalVector(U_coarse));
401
402 for (i = 0; i < n; i++)
403 {
404 u_data[i] = local_data[(int)first_index+i];
405 }
406 }
407 else
408 {
409 double *local_data;
410
411 if (my_id == 0)
412 local_data = hypre_VectorData(hypre_ParVectorLocalVector(U_coarse));
413
414 MPI_Scatterv ( local_data, info, displs, MPI_DOUBLE,
415 u_data, n, MPI_DOUBLE, 0, new_comm );
416 /*if (my_id == 0)
417 local_data = hypre_VectorData(hypre_ParVectorLocalVector(F_coarse));
418 hypre_ MPI_Scatterv ( local_data, info, displs, MPI_DOUBLE,
419 f_data, n, MPI_DOUBLE, 0, new_comm );*/
420 if (my_id == 0) hypre_TFree(displs);
421 hypre_TFree(info);
422 }
423 }
424
425 return(Solve_err_flag);
426}
427
428/* generate sub communicator, which contains only idle processors */
429
430int hypre_GenerateSubComm(MPI_Comm comm, int participate, MPI_Comm *new_comm_ptr)
431{
432 MPI_Comm new_comm;
433 MPI_Group orig_group, new_group;
434 MPI_Op hypre_MPI_MERGE;
435 int *info, *ranks, new_num_procs, my_info, my_id, num_procs;
436 int *list_len;
437
438 MPI_Comm_rank(comm,&my_id);
439
440 if (participate)
441 my_info = 1;
442 else
443 my_info = 0;
444
445 MPI_Allreduce(&my_info, &new_num_procs, 1, MPI_INT, MPI_SUM, comm);
446
447 if (new_num_procs == 0)
448 {
449 new_comm = (MPI_Comm)MPI_COMM_NULL;
450 *new_comm_ptr = new_comm;
451 return 0;
452 }
453 ranks = hypre_CTAlloc(int, new_num_procs+2);
454 if (new_num_procs == 1)
455 {
456 if (participate) my_info = my_id;
457 MPI_Allreduce(&my_info, &ranks[2], 1, MPI_INT, MPI_SUM, comm);
458 }
459 else
460 {
461 info = hypre_CTAlloc(int, new_num_procs+2);
462 list_len = hypre_CTAlloc(int, 1);
463
464 if (participate)
465 {
466 info[0] = 1;
467 info[1] = 1;
468 info[2] = my_id;
469 }
470 else
471 info[0] = 0;
472
473 list_len[0] = new_num_procs + 2;
474
475 MPI_Op_create((MPI_User_function *)hypre_merge_lists, 0, &hypre_MPI_MERGE);
476
477 MPI_Allreduce(info, ranks, list_len[0], MPI_INT, hypre_MPI_MERGE, comm);
478
479 MPI_Op_free (&hypre_MPI_MERGE);
480 hypre_TFree(list_len);
481 hypre_TFree(info);
482 }
483 MPI_Comm_size(comm,&num_procs);
484 MPI_Comm_group(comm, &orig_group);
485 MPI_Group_incl(orig_group, new_num_procs, &ranks[2], &new_group);
486 MPI_Comm_create(comm, new_group, &new_comm);
487 MPI_Group_free(&new_group);
488 MPI_Group_free(&orig_group);
489
490 hypre_TFree(ranks);
491
492 *new_comm_ptr = new_comm;
493
494 return 0;
495}
496
497
498void hypre_merge_lists (int *list1, int* list2, int *np1, MPI_Datatype *dptr)
499{
500 int i, len1, len2, indx1, indx2;
501
502 if (list1[0] == 0 || (list2[0] == 0 && list1[0] == 0))
503 {
504 return;
505 }
506 else
507 {
508 list2[0] = 1;
509 len1 = list1[1];
510 len2 = list2[1];
511 list2[1] = len1+len2;
512 if (list2[1] > *np1+2) printf("segfault in MPI User function merge_list\n");
513 indx1 = len1+1;
514 indx2 = len2+1;
515 for (i=len1+len2+1; i > 1; i--)
516 {
517 if (indx2 > 1 && indx1 > 1 && list1[indx1] > list2[indx2])
518 {
519 list2[i] = list1[indx1];
520 indx1--;
521 }
522 else if (indx2 > 1)
523 {
524 list2[i] = list2[indx2];
525 indx2--;
526 }
527 else if (indx1 > 1)
528 {
529 list2[i] = list1[indx1];
530 indx1--;
531 }
532 }
533 }
534}
535
Note: See TracBrowser for help on using the repository browser.