source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_laplace.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: 13.1 KB
RevLine 
[2aa6644]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#include "headers.h"
17
18/*--------------------------------------------------------------------------
19 * hypre_GenerateLaplacian
20 *--------------------------------------------------------------------------*/
21
22HYPRE_ParCSRMatrix
23GenerateLaplacian( MPI_Comm comm,
24 HYPRE_BigInt nx,
25 HYPRE_BigInt ny,
26 HYPRE_BigInt nz,
27 int P,
28 int Q,
29 int R,
30 int p,
31 int q,
32 int r,
33 double *value,
34 HYPRE_ParVector *rhs_ptr,
35 HYPRE_ParVector *x_ptr)
36{
37 hypre_ParCSRMatrix *A;
38 hypre_CSRMatrix *diag;
39 hypre_CSRMatrix *offd;
40 hypre_ParVector *par_rhs;
41 hypre_ParVector *par_x;
42 hypre_Vector *rhs;
43 hypre_Vector *x;
44 double *rhs_data;
45 double *x_data;
46
47 int *diag_i;
48 int *diag_j;
49 double *diag_data;
50
51 int *offd_i;
52 int *offd_j = NULL;
53 double *offd_data = NULL;
54
55 HYPRE_BigInt *global_part;
56 HYPRE_BigInt *global_part_rhs;
57 HYPRE_BigInt *global_part_x;
58 int *tmp_j = NULL;
59 int ix, iy, iz;
60 int cnt, o_cnt;
61 int local_num_rows;
62 HYPRE_BigInt *col_map_offd = NULL;
63 int row_index;
64 int i, j;
65 int ip, iq, ir;
66
67 int nx_local, ny_local, nz_local;
68 int num_cols_offd;
69 HYPRE_BigInt grid_size;
70
71 HYPRE_BigInt *nx_part;
72 HYPRE_BigInt *ny_part;
73 HYPRE_BigInt *nz_part;
74
75 int num_procs, my_id;
76 int P_busy, Q_busy, R_busy;
77
78#ifdef HYPRE_NO_GLOBAL_PARTITION
79 MPI_Comm_size(comm,&num_procs);
80 MPI_Comm_rank(comm,&my_id);
81
82 hypre_GenerateLocalPartitioning(nx,P,p,&nx_part);
83 hypre_GenerateLocalPartitioning(ny,Q,q,&ny_part);
84 hypre_GenerateLocalPartitioning(nz,R,r,&nz_part);
85
86 nx_local = (int)(nx_part[1] - nx_part[0]);
87 ny_local = (int)(ny_part[1] - ny_part[0]);
88 nz_local = (int)(nz_part[1] - nz_part[0]);
89 local_num_rows = nx_local*ny_local*nz_local;
90
91 global_part = hypre_CTAlloc(HYPRE_BigInt,2);
92 global_part_rhs = hypre_CTAlloc(HYPRE_BigInt,2);
93 global_part_x = hypre_CTAlloc(HYPRE_BigInt,2);
94
95 global_part[0] = (HYPRE_BigInt)my_id*(HYPRE_BigInt)local_num_rows;
96 global_part[1] = global_part[0]+(HYPRE_BigInt)local_num_rows;
97
98 for (i=0; i< 2; i++)
99 {
100 global_part_x[i] = global_part[i];
101 global_part_rhs[i] = global_part[i];
102 }
103 ip = p;
104 iq = q;
105 ir = r;
106#else
107 int nx_size, ny_size, nz_size;
108 MPI_Comm_size(comm,&num_procs);
109 MPI_Comm_rank(comm,&my_id);
110
111 hypre_GeneratePartitioning(nx,P,&nx_part);
112 hypre_GeneratePartitioning(ny,Q,&ny_part);
113 hypre_GeneratePartitioning(nz,R,&nz_part);
114
115 global_part = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
116 global_part_rhs = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
117 global_part_x = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
118
119 global_part[0] = 0;
120 global_part_rhs[0] = 0;
121 global_part_x[0] = 0;
122 cnt = 1;
123 for (iz = 0; iz < R; iz++)
124 {
125 nz_size = (int)(nz_part[iz+1]-nz_part[iz]);
126 for (iy = 0; iy < Q; iy++)
127 {
128 ny_size = (int)(ny_part[iy+1]-ny_part[iy]);
129 for (ix = 0; ix < P; ix++)
130 {
131 nx_size = (int)(nx_part[ix+1] - nx_part[ix]);
132 global_part[cnt] = global_part[cnt-1];
133 global_part[cnt] += (HYPRE_BigInt)(nx_size*ny_size*nz_size);
134 global_part_x[cnt] = global_part[cnt];
135 global_part_rhs[cnt] = global_part[cnt];
136 cnt++;
137 }
138 }
139 }
140
141 nx_local = (int)(nx_part[p+1] - nx_part[p]);
142 ny_local = (int)(ny_part[q+1] - ny_part[q]);
143 nz_local = (int)(nz_part[r+1] - nz_part[r]);
144 local_num_rows = nx_local*ny_local*nz_local;
145 ip = p;
146 iq = q;
147 ir = r;
148#endif
149
150 grid_size = nx*ny*nz;
151
152 diag_i = hypre_CTAlloc(int, local_num_rows+1);
153 offd_i = hypre_CTAlloc(int, local_num_rows+1);
154 rhs_data = hypre_CTAlloc(double, local_num_rows);
155 x_data = hypre_CTAlloc(double, local_num_rows);
156
157 for (i=0; i < local_num_rows; i++)
158 {
159 x_data[i] = 0.0;
160 rhs_data[i] = 1.0;
161 }
162
163 P_busy = hypre_min(nx,P);
164 Q_busy = hypre_min(ny,Q);
165 R_busy = hypre_min(nz,R);
166
167 num_cols_offd = 0;
168 if (p) num_cols_offd += ny_local*nz_local;
169 if (p < P_busy-1) num_cols_offd += ny_local*nz_local;
170 if (q) num_cols_offd += nx_local*nz_local;
171 if (q < Q_busy-1) num_cols_offd += nx_local*nz_local;
172 if (r) num_cols_offd += nx_local*ny_local;
173 if (r < R_busy-1) num_cols_offd += nx_local*ny_local;
174
175 if (!local_num_rows) num_cols_offd = 0;
176
177
178 if (num_cols_offd)
179 {
180 col_map_offd = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd);
181 }
182
183 cnt = 1;
184 o_cnt = 1;
185 diag_i[0] = 0;
186 offd_i[0] = 0;
187 for (iz = 0; iz < nz_local; iz++)
188 {
189 for (iy = 0; iy < ny_local; iy++)
190 {
191 for (ix = 0; ix < nx_local; ix++)
192 {
193 diag_i[cnt] = diag_i[cnt-1];
194 offd_i[o_cnt] = offd_i[o_cnt-1];
195 diag_i[cnt]++;
196 if (iz)
197 diag_i[cnt]++;
198 else
199 {
200 if (ir)
201 {
202 offd_i[o_cnt]++;
203 }
204 }
205 if (iy)
206 diag_i[cnt]++;
207 else
208 {
209 if (iq)
210 {
211 offd_i[o_cnt]++;
212 }
213 }
214 if (ix)
215 diag_i[cnt]++;
216 else
217 {
218 if (ip)
219 {
220 offd_i[o_cnt]++;
221 }
222 }
223 if (ix+1 < nx_local)
224 diag_i[cnt]++;
225 else
226 {
227 if (ip+1 < P)
228 {
229 offd_i[o_cnt]++;
230 }
231 }
232 if (iy+1 < ny_local)
233 diag_i[cnt]++;
234 else
235 {
236 if (iq+1 < Q)
237 {
238 offd_i[o_cnt]++;
239 }
240 }
241 if (iz+1 < nz_local)
242 diag_i[cnt]++;
243 else
244 {
245 if (ir+1 < R)
246 {
247 offd_i[o_cnt]++;
248 }
249 }
250 cnt++;
251 o_cnt++;
252 }
253 }
254 }
255
256
257 diag_j = hypre_CTAlloc(int, diag_i[local_num_rows]);
258 diag_data = hypre_CTAlloc(double, diag_i[local_num_rows]);
259
260 if (num_procs > 1)
261 {
262 offd_j = hypre_CTAlloc(int, offd_i[local_num_rows]);
263 offd_data = hypre_CTAlloc(double, offd_i[local_num_rows]);
264 }
265
266 row_index = 0;
267 cnt = 0;
268 o_cnt = 0;
269 for (iz = 0; iz < nz_local; iz++)
270 {
271 for (iy = 0; iy < ny_local; iy++)
272 {
273 for (ix = 0; ix < nx_local; ix++)
274 {
275 diag_j[cnt] = row_index;
276 diag_data[cnt++] = value[0];
277 if (iz)
278 {
279 diag_j[cnt] = row_index-nx_local*ny_local;
280 diag_data[cnt++] = value[3];
281 }
282 else
283 {
284 if (ir)
285 {
286 offd_j[o_cnt] = o_cnt;
287 hypre_map(ix,iy,iz-1,p,q,r-1,P,Q,R,
288 nx_part,ny_part,nz_part,global_part,
289 &col_map_offd[o_cnt]);
290 offd_data[o_cnt++] = value[3];
291 }
292 }
293 if (iy)
294 {
295 diag_j[cnt] = row_index-nx_local;
296 diag_data[cnt++] = value[2];
297 }
298 else
299 {
300 if (iq)
301 {
302 offd_j[o_cnt] = o_cnt;
303 hypre_map(ix,iy-1,iz,p,q-1,r,P,Q,R,
304 nx_part,ny_part,nz_part,global_part,
305 &col_map_offd[o_cnt]);
306 offd_data[o_cnt++] = value[2];
307 }
308 }
309 if (ix)
310 {
311 diag_j[cnt] = row_index-1;
312 diag_data[cnt++] = value[1];
313 }
314 else
315 {
316 if (ip)
317 {
318 offd_j[o_cnt] = o_cnt;
319 hypre_map(ix-1,iy,iz,p-1,q,r,P,Q,R,
320 nx_part,ny_part,nz_part,global_part,
321 &col_map_offd[o_cnt]);
322 offd_data[o_cnt++] = value[1];
323 }
324 }
325 if (ix+1 < nx_local)
326 {
327 diag_j[cnt] = row_index+1;
328 diag_data[cnt++] = value[1];
329 }
330 else
331 {
332 if (ip+1 < P)
333 {
334 offd_j[o_cnt] = o_cnt;
335 hypre_map(ix+1,iy,iz,p+1,q,r,P,Q,R,
336 nx_part,ny_part,nz_part,global_part,
337 &col_map_offd[o_cnt]);
338 offd_data[o_cnt++] = value[1];
339 }
340 }
341 if (iy+1 < ny_local)
342 {
343 diag_j[cnt] = row_index+nx_local;
344 diag_data[cnt++] = value[2];
345 }
346 else
347 {
348 if (iq+1 < Q)
349 {
350 offd_j[o_cnt] = o_cnt;
351 hypre_map(ix,iy+1,iz,p,q+1,r,P,Q,R,
352 nx_part,ny_part,nz_part,global_part,
353 &col_map_offd[o_cnt]);
354 offd_data[o_cnt++] = value[2];
355 }
356 }
357 if (iz+1 < nz_local)
358 {
359 diag_j[cnt] = row_index+nx_local*ny_local;
360 diag_data[cnt++] = value[3];
361 }
362 else
363 {
364 if (ir+1 < R)
365 {
366 offd_j[o_cnt] = o_cnt;
367 hypre_map(ix,iy,iz+1,p,q,r+1,P,Q,R,
368 nx_part,ny_part,nz_part,global_part,
369 &col_map_offd[o_cnt]);
370 offd_data[o_cnt++] = value[3];
371 }
372 }
373 row_index++;
374 }
375 }
376 }
377
378
379 if (num_cols_offd) tmp_j = hypre_CTAlloc(int, num_cols_offd);
380 for (i=0; i < num_cols_offd; i++)
381 tmp_j[i] = offd_j[i];
382
383
384 if (num_procs > 1)
385 {
386 hypre_BigQsortbi(col_map_offd, tmp_j, 0, num_cols_offd-1);
387
388 for (i=0; i < num_cols_offd; i++)
389 for (j=0; j < num_cols_offd; j++)
390 if (offd_j[i] == tmp_j[j])
391 {
392 offd_j[i] = j;
393 break;
394 }
395 }
396
397
398 hypre_TFree(tmp_j);
399
400 A = hypre_ParCSRMatrixCreate(comm, grid_size, grid_size,
401 global_part, global_part, num_cols_offd,
402 diag_i[local_num_rows],
403 offd_i[local_num_rows]);
404
405 par_rhs = hypre_ParVectorCreate(comm, grid_size, global_part_rhs);
406 rhs = hypre_ParVectorLocalVector(par_rhs);
407 hypre_VectorData(rhs) = rhs_data;
408
409 par_x = hypre_ParVectorCreate(comm, grid_size, global_part_x);
410 x = hypre_ParVectorLocalVector(par_x);
411 hypre_VectorData(x) = x_data;
412
413 hypre_ParCSRMatrixColMapOffd(A) = col_map_offd;
414
415 diag = hypre_ParCSRMatrixDiag(A);
416 hypre_CSRMatrixI(diag) = diag_i;
417 hypre_CSRMatrixJ(diag) = diag_j;
418 hypre_CSRMatrixData(diag) = diag_data;
419
420 offd = hypre_ParCSRMatrixOffd(A);
421 hypre_CSRMatrixI(offd) = offd_i;
422 if (num_cols_offd)
423 {
424 hypre_CSRMatrixJ(offd) = offd_j;
425 hypre_CSRMatrixData(offd) = offd_data;
426 }
427
428 hypre_TFree(nx_part);
429 hypre_TFree(ny_part);
430 hypre_TFree(nz_part);
431
432 *rhs_ptr = (HYPRE_ParVector) par_rhs;
433 *x_ptr = (HYPRE_ParVector) par_x;
434
435 return (HYPRE_ParCSRMatrix) A;
436}
437
438/*--------------------------------------------------------------------------
439 *--------------------------------------------------------------------------*/
440
441int
442hypre_map( int ix,
443 int iy,
444 int iz,
445 int p,
446 int q,
447 int r,
448 int P,
449 int Q,
450 int R,
451 HYPRE_BigInt *nx_part,
452 HYPRE_BigInt *ny_part,
453 HYPRE_BigInt *nz_part,
454 HYPRE_BigInt *global_part,
455 HYPRE_BigInt *value_ptr)
456{
457 int nx_local;
458 int ny_local;
459 int nz_local;
460 HYPRE_BigInt global_index;
461 int proc_num;
462
463#ifdef HYPRE_NO_GLOBAL_PARTITION
464 proc_num = r*P*Q + q*P + p;
465 nx_local = (int)(nx_part[1] - nx_part[0]);
466 ny_local = (int)(ny_part[1] - ny_part[0]);
467 nz_local = (int)(nz_part[1] - nz_part[0]);
468 if (ix < 0) ix += nx_local;
469 if (ix >= nx_local) ix -= nx_local;
470 if (iy < 0) iy += ny_local;
471 if (iy >= ny_local) iy -= ny_local;
472 if (iz < 0) iz += nz_local;
473 if (iz >= nz_local) iz -= nz_local;
474 global_index = (HYPRE_BigInt)proc_num
475 *(HYPRE_BigInt)(nx_local*ny_local*nz_local)
476 + (HYPRE_BigInt)((iz*ny_local+iy)*nx_local + ix);
477#else
478 proc_num = r*P*Q + q*P + p;
479 nx_local = (int)(nx_part[p+1] - nx_part[p]);
480 ny_local = (int)(ny_part[q+1] - ny_part[q]);
481 nz_local = (int)(nz_part[r+1] - nz_part[r]);
482 if (ix < 0) ix += nx_local;
483 if (ix >= nx_local) ix = 0;
484 if (iy < 0) iy += ny_local;
485 if (iy >= ny_local) iy = 0;
486 if (iz < 0) iz += nz_local;
487 if (iz >= nz_local) iz = 0;
488 global_index = global_part[proc_num]
489 + (HYPRE_BigInt)((iz*ny_local+iy)*nx_local + ix);
490#endif
491 *value_ptr = global_index;
492
493 return 0;
494}
Note: See TracBrowser for help on using the repository browser.