source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_nodal_systems.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: 31.2 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 *
17 *****************************************************************************/
18
19/* following should be in a header file */
20
21
22#include "headers.h"
23
24
25
26/*==========================================================================*/
27/*==========================================================================*/
28/**
29 Generates nodal norm matrix for use with nodal systems version
30
31 {\bf Input files:}
32 headers.h
33
34 @return Error code.
35
36 @param A [IN]
37 coefficient matrix
38 @param AN_ptr [OUT]
39 nodal norm matrix
40
41 @see */
42/*--------------------------------------------------------------------------*/
43
44int
45hypre_BoomerAMGCreateNodalA(hypre_ParCSRMatrix *A,
46 int num_functions,
47 int *dof_func,
48 int option,
49 hypre_ParCSRMatrix **AN_ptr)
50{
51 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
52 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
53 int *A_diag_i = hypre_CSRMatrixI(A_diag);
54 double *A_diag_data = hypre_CSRMatrixData(A_diag);
55
56
57 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
58 int *A_offd_i = hypre_CSRMatrixI(A_offd);
59 double *A_offd_data = hypre_CSRMatrixData(A_offd);
60 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
61 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
62
63 HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(A);
64 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
65 int num_variables = hypre_CSRMatrixNumRows(A_diag);
66 int num_nonzeros_offd = 0;
67 int num_cols_offd = 0;
68
69 hypre_ParCSRMatrix *AN;
70 hypre_CSRMatrix *AN_diag;
71 int *AN_diag_i;
72 int *AN_diag_j;
73 double *AN_diag_data;
74 hypre_CSRMatrix *AN_offd;
75 int *AN_offd_i;
76 int *AN_offd_j = NULL;
77 double *AN_offd_data = NULL;
78 HYPRE_BigInt *col_map_offd_AN = NULL;
79 HYPRE_BigInt *new_col_map_offd = NULL;
80 HYPRE_BigInt *row_starts_AN;
81 int AN_num_nonzeros_diag = 0;
82 int AN_num_nonzeros_offd = 0;
83 int num_cols_offd_AN;
84 int new_num_cols_offd;
85
86 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
87 int num_sends;
88 int num_recvs;
89 int *send_procs;
90 int *send_map_starts;
91 int *send_map_elmts;
92 int *new_send_map_elmts;
93 int *recv_procs;
94 int *recv_vec_starts;
95
96 hypre_ParCSRCommPkg *comm_pkg_AN;
97 int *send_procs_AN;
98 int *send_map_starts_AN;
99 int *send_map_elmts_AN;
100 int *recv_procs_AN;
101 int *recv_vec_starts_AN;
102
103 int i, j, k, k_map;
104
105 int ierr = 0;
106
107 int index, row;
108 int start_index;
109 int num_procs;
110 int mode, node, cnt;
111 HYPRE_BigInt big_node;
112 int new_send_elmts_size;
113
114 HYPRE_BigInt global_num_nodes;
115 int num_nodes;
116 int num_fun2;
117 HYPRE_BigInt *big_map_to_node = NULL;
118 int *map_to_node = NULL;
119 int *map_to_map;
120 int *counter;
121
122 MPI_Comm_size(comm,&num_procs);
123
124 if (!comm_pkg)
125 {
126#ifdef HYPRE_NO_GLOBAL_PARTITION
127 hypre_NewCommPkgCreate(A);
128#else
129 hypre_MatvecCommPkgCreate(A);
130#endif
131 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
132 }
133
134 mode = fabs(option);
135
136 comm_pkg_AN = NULL;
137 col_map_offd_AN = NULL;
138
139#ifdef HYPRE_NO_GLOBAL_PARTITION
140 row_starts_AN = hypre_CTAlloc(HYPRE_BigInt, 2);
141
142 for (i=0; i < 2; i++)
143 {
144 row_starts_AN[i] = row_starts[i]/(HYPRE_BigInt)num_functions;
145 if (row_starts_AN[i]*(HYPRE_BigInt)num_functions < row_starts[i])
146 {
147 printf("nodes not properly aligned or incomplete info!\n");
148 return (87);
149 }
150 }
151
152 global_num_nodes = hypre_ParCSRMatrixGlobalNumRows(A)/(HYPRE_BigInt)num_functions;
153
154
155#else
156 row_starts_AN = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
157
158 for (i=0; i < num_procs+1; i++)
159 {
160 row_starts_AN[i] = row_starts[i]/(HYPRE_BigInt)num_functions;
161 if (row_starts_AN[i]*(HYPRE_BigInt)num_functions < row_starts[i])
162 {
163 printf("nodes not properly aligned or incomplete info!\n");
164 return (87);
165 }
166 }
167
168 global_num_nodes = row_starts_AN[num_procs];
169
170#endif
171
172
173 num_nodes = num_variables/num_functions;
174 num_fun2 = num_functions*num_functions;
175
176 map_to_node = hypre_CTAlloc(int, num_variables);
177 AN_diag_i = hypre_CTAlloc(int, num_nodes+1);
178 counter = hypre_CTAlloc(int, num_nodes);
179 for (i=0; i < num_variables; i++)
180 map_to_node[i] = i/num_functions;
181 for (i=0; i < num_nodes; i++)
182 counter[i] = -1;
183
184 AN_num_nonzeros_diag = 0;
185 row = 0;
186 for (i=0; i < num_nodes; i++)
187 {
188 AN_diag_i[i] = AN_num_nonzeros_diag;
189 for (j=0; j < num_functions; j++)
190 {
191 for (k=A_diag_i[row]; k < A_diag_i[row+1]; k++)
192 {
193 k_map = map_to_node[A_diag_j[k]];
194 if (counter[k_map] < i)
195 {
196 counter[k_map] = i;
197 AN_num_nonzeros_diag++;
198 }
199 }
200 row++;
201 }
202 }
203 AN_diag_i[num_nodes] = AN_num_nonzeros_diag;
204
205 AN_diag_j = hypre_CTAlloc(int, AN_num_nonzeros_diag);
206 AN_diag_data = hypre_CTAlloc(double, AN_num_nonzeros_diag);
207
208 AN_diag = hypre_CSRMatrixCreate(num_nodes,num_nodes,AN_num_nonzeros_diag);
209 hypre_CSRMatrixI(AN_diag) = AN_diag_i;
210 hypre_CSRMatrixJ(AN_diag) = AN_diag_j;
211 hypre_CSRMatrixData(AN_diag) = AN_diag_data;
212
213 for (i=0; i < num_nodes; i++)
214 counter[i] = -1;
215 index = 0;
216 start_index = 0;
217 row = 0;
218
219 switch (mode)
220 {
221 case 7: /* frobenius norm with signs*/
222 {
223 for (i=0; i < num_nodes; i++)
224 {
225 for (j=0; j < num_functions; j++)
226 {
227 for (k=A_diag_i[row]; k < A_diag_i[row+1]; k++)
228 {
229 k_map = map_to_node[A_diag_j[k]];
230 if (counter[k_map] < start_index)
231 {
232 counter[k_map] = index;
233 AN_diag_j[index] = k_map;
234 AN_diag_data[index] = A_diag_data[k]*A_diag_data[k];
235 index++;
236 }
237 else
238 {
239 AN_diag_data[counter[k_map]] +=
240 A_diag_data[k]*A_diag_data[k];
241 }
242 }
243 row++;
244 }
245 start_index = index;
246 }
247 for (i=0; i < AN_num_nonzeros_diag; i++)
248 AN_diag_data[i] = sqrt(AN_diag_data[i]);
249
250 /* temp for testing - make all diagonal entries negative */
251 /* the diagonal is the first element listed in each row -
252 this is the same as serial code */
253
254 for (i=0; i < num_nodes; i++)
255 {
256 index = AN_diag_i[i];
257 AN_diag_data[index] = - AN_diag_data[index];
258 }
259
260 }
261 break;
262
263 case 1: /* frobenius norm */
264 {
265 for (i=0; i < num_nodes; i++)
266 {
267 for (j=0; j < num_functions; j++)
268 {
269 for (k=A_diag_i[row]; k < A_diag_i[row+1]; k++)
270 {
271 k_map = map_to_node[A_diag_j[k]];
272 if (counter[k_map] < start_index)
273 {
274 counter[k_map] = index;
275 AN_diag_j[index] = k_map;
276 AN_diag_data[index] = A_diag_data[k]*A_diag_data[k];
277 index++;
278 }
279 else
280 {
281 AN_diag_data[counter[k_map]] +=
282 A_diag_data[k]*A_diag_data[k];
283 }
284 }
285 row++;
286 }
287 start_index = index;
288 }
289 for (i=0; i < AN_num_nonzeros_diag; i++)
290 AN_diag_data[i] = sqrt(AN_diag_data[i]);
291
292#if 0
293 /* temp for testing - make all diagonal entries negative */
294 /* the diagonal is the first element listed in each row -
295 this is the same as serial code */
296
297 for (i=0; i < num_nodes; i++)
298 {
299 index = AN_diag_i[i];
300 AN_diag_data[index] = - AN_diag_data[index];
301 }
302#endif
303
304 }
305 break;
306
307 case 2: /* sum of abs. value of all elements in each block */
308 {
309 for (i=0; i < num_nodes; i++)
310 {
311 for (j=0; j < num_functions; j++)
312 {
313 for (k=A_diag_i[row]; k < A_diag_i[row+1]; k++)
314 {
315 k_map = map_to_node[A_diag_j[k]];
316 if (counter[k_map] < start_index)
317 {
318 counter[k_map] = index;
319 AN_diag_j[index] = k_map;
320 AN_diag_data[index] = fabs(A_diag_data[k]);
321 index++;
322 }
323 else
324 {
325 AN_diag_data[counter[k_map]] += fabs(A_diag_data[k]);
326 }
327 }
328 row++;
329 }
330 start_index = index;
331 }
332 for (i=0; i < AN_num_nonzeros_diag; i++)
333 AN_diag_data[i] /= num_fun2;
334 }
335 break;
336
337 case 3: /* largest element of each block (sets true value - not abs. value) */
338 {
339
340 for (i=0; i < num_nodes; i++)
341 {
342 for (j=0; j < num_functions; j++)
343 {
344 for (k=A_diag_i[row]; k < A_diag_i[row+1]; k++)
345 {
346 k_map = map_to_node[A_diag_j[k]];
347 if (counter[k_map] < start_index)
348 {
349 counter[k_map] = index;
350 AN_diag_j[index] = k_map;
351 AN_diag_data[index] = A_diag_data[k];
352 index++;
353 }
354 else
355 {
356 if (fabs(A_diag_data[k]) >
357 fabs(AN_diag_data[counter[k_map]]))
358 AN_diag_data[counter[k_map]] = A_diag_data[k];
359 }
360 }
361 row++;
362 }
363 start_index = index;
364 }
365 }
366 break;
367 }
368
369 num_nonzeros_offd = A_offd_i[num_variables];
370 AN_offd_i = hypre_CTAlloc(int, num_nodes+1);
371
372 num_cols_offd_AN = 0;
373
374 if (comm_pkg)
375 {
376 comm_pkg_AN = hypre_CTAlloc(hypre_ParCSRCommPkg,1);
377 hypre_ParCSRCommPkgComm(comm_pkg_AN) = comm;
378 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
379 hypre_ParCSRCommPkgNumSends(comm_pkg_AN) = num_sends;
380 num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
381 hypre_ParCSRCommPkgNumRecvs(comm_pkg_AN) = num_recvs;
382 send_procs = hypre_ParCSRCommPkgSendProcs(comm_pkg);
383 send_map_starts = hypre_ParCSRCommPkgSendMapStarts(comm_pkg);
384 send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
385 recv_procs = hypre_ParCSRCommPkgRecvProcs(comm_pkg);
386 recv_vec_starts = hypre_ParCSRCommPkgRecvVecStarts(comm_pkg);
387 send_procs_AN = NULL;
388 send_map_elmts_AN = NULL;
389 if (num_sends)
390 {
391 send_procs_AN = hypre_CTAlloc(int,num_sends);
392 send_map_elmts_AN = hypre_CTAlloc(int,send_map_starts[num_sends]);
393 }
394 send_map_starts_AN = hypre_CTAlloc(int,num_sends+1);
395 recv_vec_starts_AN = hypre_CTAlloc(int,num_recvs+1);
396 recv_procs_AN = NULL;
397 if (num_recvs) recv_procs_AN = hypre_CTAlloc(int,num_recvs);
398 for (i=0; i < num_sends; i++)
399 send_procs_AN[i] = send_procs[i];
400 for (i=0; i < num_recvs; i++)
401 recv_procs_AN[i] = recv_procs[i];
402
403 send_map_starts_AN[0] = 0;
404 cnt = 0;
405 for (i=0; i < num_sends; i++)
406 {
407 k_map = send_map_starts[i];
408 if (send_map_starts[i+1]-k_map)
409 send_map_elmts_AN[cnt++] = send_map_elmts[k_map]/num_functions;
410 for (j=send_map_starts[i]+1; j < send_map_starts[i+1]; j++)
411 {
412 node = send_map_elmts[j]/num_functions;
413 if (node > send_map_elmts_AN[cnt-1])
414 send_map_elmts_AN[cnt++] = node;
415 }
416 send_map_starts_AN[i+1] = cnt;
417 }
418 hypre_ParCSRCommPkgSendProcs(comm_pkg_AN) = send_procs_AN;
419 hypre_ParCSRCommPkgSendMapStarts(comm_pkg_AN) = send_map_starts_AN;
420 hypre_ParCSRCommPkgSendMapElmts(comm_pkg_AN) = send_map_elmts_AN;
421 hypre_ParCSRCommPkgRecvProcs(comm_pkg_AN) = recv_procs_AN;
422 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_AN) = recv_vec_starts_AN;
423 }
424
425 num_cols_offd = hypre_CSRMatrixNumCols(A_offd);
426 hypre_TFree(map_to_node);
427 if (num_cols_offd)
428 {
429 if (num_cols_offd > num_variables)
430 {
431 big_map_to_node = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd);
432 }
433
434 num_cols_offd_AN = 1;
435 big_map_to_node[0] = col_map_offd[0]/(HYPRE_BigInt)num_functions;
436 for (i=1; i < num_cols_offd; i++)
437 {
438 big_map_to_node[i] = col_map_offd[i]/(HYPRE_BigInt)num_functions;
439 if (big_map_to_node[i] > big_map_to_node[i-1]) num_cols_offd_AN++;
440 }
441
442 if (num_cols_offd_AN > num_nodes)
443 {
444 hypre_TFree(counter);
445 counter = hypre_CTAlloc(int,num_cols_offd_AN);
446 }
447
448 map_to_map = NULL;
449 col_map_offd_AN = NULL;
450 map_to_map = hypre_CTAlloc(int, num_cols_offd);
451 col_map_offd_AN = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd_AN);
452 col_map_offd_AN[0] = big_map_to_node[0];
453 recv_vec_starts_AN[0] = 0;
454 cnt = 1;
455 for (i=0; i < num_recvs; i++)
456 {
457 for (j=recv_vec_starts[i]; j < recv_vec_starts[i+1]; j++)
458 {
459 big_node = big_map_to_node[j];
460 if (big_node > col_map_offd_AN[cnt-1])
461 {
462 col_map_offd_AN[cnt++] = big_node;
463 }
464 map_to_map[j] = cnt-1;
465 }
466 recv_vec_starts_AN[i+1] = cnt;
467 }
468
469 for (i=0; i < num_cols_offd_AN; i++)
470 counter[i] = -1;
471
472 AN_num_nonzeros_offd = 0;
473 row = 0;
474 for (i=0; i < num_nodes; i++)
475 {
476 AN_offd_i[i] = AN_num_nonzeros_offd;
477 for (j=0; j < num_functions; j++)
478 {
479 for (k=A_offd_i[row]; k < A_offd_i[row+1]; k++)
480 {
481 k_map = map_to_map[A_offd_j[k]];
482 if (counter[k_map] < i)
483 {
484 counter[k_map] = i;
485 AN_num_nonzeros_offd++;
486 }
487 }
488 row++;
489 }
490 }
491 AN_offd_i[num_nodes] = AN_num_nonzeros_offd;
492 }
493
494
495 AN_offd = hypre_CSRMatrixCreate(num_nodes,num_cols_offd_AN,
496 AN_num_nonzeros_offd);
497 hypre_CSRMatrixI(AN_offd) = AN_offd_i;
498 if (AN_num_nonzeros_offd)
499 {
500 AN_offd_j = hypre_CTAlloc(int, AN_num_nonzeros_offd);
501 AN_offd_data = hypre_CTAlloc(double, AN_num_nonzeros_offd);
502 hypre_CSRMatrixJ(AN_offd) = AN_offd_j;
503 hypre_CSRMatrixData(AN_offd) = AN_offd_data;
504
505 for (i=0; i < num_cols_offd_AN; i++)
506 counter[i] = -1;
507 index = 0;
508 row = 0;
509 AN_offd_i[0] = 0;
510 start_index = 0;
511 switch (mode)
512 {
513 case 1: /* frobenius norm */
514 {
515 for (i=0; i < num_nodes; i++)
516 {
517 for (j=0; j < num_functions; j++)
518 {
519 for (k=A_offd_i[row]; k < A_offd_i[row+1]; k++)
520 {
521 k_map = map_to_map[A_offd_j[k]];
522 if (counter[k_map] < start_index)
523 {
524 counter[k_map] = index;
525 AN_offd_j[index] = k_map;
526 AN_offd_data[index] = A_offd_data[k]*A_offd_data[k];
527 index++;
528 }
529 else
530 {
531 AN_offd_data[counter[k_map]] +=
532 A_offd_data[k]*A_offd_data[k];
533 }
534 }
535 row++;
536 }
537 start_index = index;
538 }
539 for (i=0; i < AN_num_nonzeros_offd; i++)
540 AN_offd_data[i] = sqrt(AN_offd_data[i]);
541 }
542 break;
543
544 case 2: /* sum of abs. value of all elements in block */
545 {
546 for (i=0; i < num_nodes; i++)
547 {
548 for (j=0; j < num_functions; j++)
549 {
550 for (k=A_offd_i[row]; k < A_offd_i[row+1]; k++)
551 {
552 k_map = map_to_map[A_offd_j[k]];
553 if (counter[k_map] < start_index)
554 {
555 counter[k_map] = index;
556 AN_offd_j[index] = k_map;
557 AN_offd_data[index] = fabs(A_offd_data[k]);
558 index++;
559 }
560 else
561 {
562 AN_offd_data[counter[k_map]] += fabs(A_offd_data[k]);
563 }
564 }
565 row++;
566 }
567 start_index = index;
568 }
569 for (i=0; i < AN_num_nonzeros_offd; i++)
570 AN_offd_data[i] /= num_fun2;
571 }
572 break;
573
574 case 3: /* largest element in each block (not abs. value ) */
575 {
576 for (i=0; i < num_nodes; i++)
577 {
578 for (j=0; j < num_functions; j++)
579 {
580 for (k=A_offd_i[row]; k < A_offd_i[row+1]; k++)
581 {
582 k_map = map_to_map[A_offd_j[k]];
583 if (counter[k_map] < start_index)
584 {
585 counter[k_map] = index;
586 AN_offd_j[index] = k_map;
587 AN_offd_data[index] = A_offd_data[k];
588 index++;
589 }
590 else
591 {
592 if (fabs(A_offd_data[k]) >
593 fabs(AN_offd_data[counter[k_map]]))
594 AN_offd_data[counter[k_map]] = A_offd_data[k];
595 }
596 }
597 row++;
598 }
599 start_index = index;
600 }
601 }
602 break;
603 }
604 hypre_TFree(map_to_map);
605 }
606
607
608 AN = hypre_ParCSRMatrixCreate(comm, global_num_nodes, global_num_nodes,
609 row_starts_AN, row_starts_AN, num_cols_offd_AN,
610 AN_num_nonzeros_diag, AN_num_nonzeros_offd);
611
612 /* we already created the diag and offd matrices - so we don't need the ones
613 created above */
614 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(AN));
615 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(AN));
616 hypre_ParCSRMatrixDiag(AN) = AN_diag;
617 hypre_ParCSRMatrixOffd(AN) = AN_offd;
618
619
620 hypre_ParCSRMatrixColMapOffd(AN) = col_map_offd_AN;
621 hypre_ParCSRMatrixCommPkg(AN) = comm_pkg_AN;
622
623 new_num_cols_offd = num_functions*num_cols_offd_AN;
624
625 if (new_num_cols_offd > num_cols_offd)
626 {
627 new_col_map_offd = hypre_CTAlloc(HYPRE_BigInt, new_num_cols_offd);
628 cnt = 0;
629 for (i=0; i < num_cols_offd_AN; i++)
630 {
631 for (j=0; j < num_functions; j++)
632 {
633 new_col_map_offd[cnt++] = (HYPRE_BigInt)num_functions*
634 col_map_offd_AN[i]+(HYPRE_BigInt)j;
635 }
636 }
637 cnt = 0;
638 for (i=0; i < num_cols_offd; i++)
639 {
640 while (col_map_offd[i] > new_col_map_offd[cnt])
641 cnt++;
642 col_map_offd[i] = cnt++;
643 }
644 for (i=0; i < num_recvs+1; i++)
645 {
646 recv_vec_starts[i] = num_functions*recv_vec_starts_AN[i];
647 }
648
649 for (i=0; i < num_nonzeros_offd; i++)
650 {
651 j = A_offd_j[i];
652 A_offd_j[i] = col_map_offd[j];
653 }
654 hypre_ParCSRMatrixColMapOffd(A) = new_col_map_offd;
655 hypre_CSRMatrixNumCols(A_offd) = new_num_cols_offd;
656 hypre_TFree(col_map_offd);
657 }
658
659 hypre_TFree(big_map_to_node);
660 new_send_elmts_size = send_map_starts_AN[num_sends]*num_functions;
661
662 if (new_send_elmts_size > send_map_starts[num_sends])
663 {
664 new_send_map_elmts = hypre_CTAlloc(int,new_send_elmts_size);
665 cnt = 0;
666 send_map_starts[0] = 0;
667 for (i=0; i < num_sends; i++)
668 {
669 send_map_starts[i+1] = send_map_starts_AN[i+1]*num_functions;
670 for (j=send_map_starts_AN[i]; j < send_map_starts_AN[i+1]; j++)
671 {
672 for (k=0; k < num_functions; k++)
673 new_send_map_elmts[cnt++] = send_map_elmts_AN[j]*num_functions+k;
674 }
675 }
676 hypre_TFree(send_map_elmts);
677 hypre_ParCSRCommPkgSendMapElmts(comm_pkg) = new_send_map_elmts;
678 }
679
680 *AN_ptr = AN;
681
682 hypre_TFree(counter);
683
684 return (ierr);
685}
686
687
688/* This creates a scalar version of the CF_marker, dof_array and strength matrix (SN) */
689
690int
691hypre_BoomerAMGCreateScalarCFS(hypre_ParCSRMatrix *SN,
692 int *CFN_marker,
693 int *col_offd_SN_to_AN,
694 int num_functions,
695 int nodal,
696 int data,
697 int **dof_func_ptr,
698 int **CF_marker_ptr,
699 int **col_offd_S_to_A_ptr,
700 hypre_ParCSRMatrix **S_ptr)
701{
702 MPI_Comm comm = hypre_ParCSRMatrixComm(SN);
703 hypre_ParCSRMatrix *S;
704 hypre_CSRMatrix *S_diag;
705 int *S_diag_i;
706 int *S_diag_j;
707 double *S_diag_data;
708 hypre_CSRMatrix *S_offd;
709 int *S_offd_i;
710 int *S_offd_j;
711 double *S_offd_data;
712 HYPRE_BigInt *row_starts_S;
713 HYPRE_BigInt *col_starts_S;
714 HYPRE_BigInt *row_starts_SN = hypre_ParCSRMatrixRowStarts(SN);
715 HYPRE_BigInt *col_starts_SN = hypre_ParCSRMatrixColStarts(SN);
716 hypre_CSRMatrix *SN_diag = hypre_ParCSRMatrixDiag(SN);
717 int *SN_diag_i = hypre_CSRMatrixI(SN_diag);
718 int *SN_diag_j = hypre_CSRMatrixJ(SN_diag);
719 double *SN_diag_data;
720 hypre_CSRMatrix *SN_offd = hypre_ParCSRMatrixOffd(SN);
721 int *SN_offd_i = hypre_CSRMatrixI(SN_offd);
722 int *SN_offd_j = hypre_CSRMatrixJ(SN_offd);
723 double *SN_offd_data;
724 int *CF_marker;
725 HYPRE_BigInt *col_map_offd_SN = hypre_ParCSRMatrixColMapOffd(SN);
726 HYPRE_BigInt *col_map_offd_S;
727 int *dof_func;
728 int num_nodes = hypre_CSRMatrixNumRows(SN_diag);
729 int num_variables;
730 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(SN);
731 int num_sends;
732 int num_recvs;
733 int *send_procs;
734 int *send_map_starts;
735 int *send_map_elmts;
736 int *recv_procs;
737 int *recv_vec_starts;
738 hypre_ParCSRCommPkg *comm_pkg_S;
739 int *send_procs_S;
740 int *send_map_starts_S;
741 int *send_map_elmts_S;
742 int *recv_procs_S;
743 int *recv_vec_starts_S;
744 int *col_offd_S_to_A = NULL;
745
746 int num_coarse_nodes;
747 int i,j,k,k1,jj,cnt;
748 int row, start, end;
749 int num_procs;
750 int num_cols_offd_SN = hypre_CSRMatrixNumCols(SN_offd);
751 int num_cols_offd_S;
752 int SN_num_nonzeros_diag;
753 int SN_num_nonzeros_offd;
754 int S_num_nonzeros_diag;
755 int S_num_nonzeros_offd;
756 HYPRE_BigInt global_num_vars;
757 HYPRE_BigInt global_num_cols;
758 HYPRE_BigInt global_num_nodes;
759 int ierr = 0;
760
761 MPI_Comm_size(comm, &num_procs);
762
763 num_variables = num_functions*num_nodes;
764 CF_marker = hypre_CTAlloc(int, num_variables);
765
766 if (nodal < 0)
767 {
768 cnt = 0;
769 num_coarse_nodes = 0;
770 for (i=0; i < num_nodes; i++)
771 {
772 if (CFN_marker[i] == 1) num_coarse_nodes++;
773 for (j=0; j < num_functions; j++)
774 CF_marker[cnt++] = CFN_marker[i];
775 }
776
777 dof_func = hypre_CTAlloc(int,num_coarse_nodes*num_functions);
778 cnt = 0;
779 for (i=0; i < num_nodes; i++)
780 {
781 if (CFN_marker[i] == 1)
782 {
783 for (k=0; k < num_functions; k++)
784 dof_func[cnt++] = k;
785 }
786 }
787 *dof_func_ptr = dof_func;
788 }
789 else
790 {
791 cnt = 0;
792 for (i=0; i < num_nodes; i++)
793 for (j=0; j < num_functions; j++)
794 CF_marker[cnt++] = CFN_marker[i];
795 }
796
797 *CF_marker_ptr = CF_marker;
798
799
800#ifdef HYPRE_NO_GLOBAL_PARTITION
801 row_starts_S = hypre_CTAlloc(HYPRE_BigInt,2);
802 for (i=0; i < 2; i++)
803 row_starts_S[i] = num_functions*row_starts_SN[i];
804
805 if (row_starts_SN != col_starts_SN)
806 {
807 col_starts_S = hypre_CTAlloc(HYPRE_BigInt,2);
808 for (i=0; i < 2; i++)
809 col_starts_S[i] = num_functions*col_starts_SN[i];
810 }
811 else
812 {
813 col_starts_S = row_starts_S;
814 }
815#else
816 row_starts_S = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
817 for (i=0; i < num_procs+1; i++)
818 row_starts_S[i] = num_functions*row_starts_SN[i];
819
820 if (row_starts_SN != col_starts_SN)
821 {
822 col_starts_S = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
823 for (i=0; i < num_procs+1; i++)
824 col_starts_S[i] = num_functions*col_starts_SN[i];
825 }
826 else
827 {
828 col_starts_S = row_starts_S;
829 }
830#endif
831
832
833 SN_num_nonzeros_diag = SN_diag_i[num_nodes];
834 SN_num_nonzeros_offd = SN_offd_i[num_nodes];
835
836 global_num_nodes = hypre_ParCSRMatrixGlobalNumRows(SN);
837 global_num_cols = hypre_ParCSRMatrixGlobalNumCols(SN)*num_functions;
838
839 global_num_vars = global_num_nodes*num_functions;
840 S_num_nonzeros_diag = num_functions*SN_num_nonzeros_diag;
841 S_num_nonzeros_offd = num_functions*SN_num_nonzeros_offd;
842 num_cols_offd_S = num_functions*num_cols_offd_SN;
843 S = hypre_ParCSRMatrixCreate(comm, global_num_vars, global_num_cols,
844 row_starts_S, col_starts_S, num_cols_offd_S,
845 S_num_nonzeros_diag, S_num_nonzeros_offd);
846
847 S_diag = hypre_ParCSRMatrixDiag(S);
848 S_offd = hypre_ParCSRMatrixOffd(S);
849 S_diag_i = hypre_CTAlloc(int, num_variables+1);
850 S_offd_i = hypre_CTAlloc(int, num_variables+1);
851 S_diag_j = hypre_CTAlloc(int, S_num_nonzeros_diag);
852 hypre_CSRMatrixI(S_diag) = S_diag_i;
853 hypre_CSRMatrixJ(S_diag) = S_diag_j;
854 if (data)
855 {
856 SN_diag_data = hypre_CSRMatrixData(SN_diag);
857 S_diag_data = hypre_CTAlloc(double, S_num_nonzeros_diag);
858 hypre_CSRMatrixData(S_diag) = S_diag_data;
859 if (num_cols_offd_S)
860 {
861 SN_offd_data = hypre_CSRMatrixData(SN_offd);
862 S_offd_data = hypre_CTAlloc(double, S_num_nonzeros_offd);
863 hypre_CSRMatrixData(S_offd) = S_offd_data;
864 }
865
866 }
867 hypre_CSRMatrixI(S_offd) = S_offd_i;
868
869 if (comm_pkg)
870 {
871 comm_pkg_S = hypre_CTAlloc(hypre_ParCSRCommPkg,1);
872 hypre_ParCSRCommPkgComm(comm_pkg_S) = comm;
873 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
874 hypre_ParCSRCommPkgNumSends(comm_pkg_S) = num_sends;
875 num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
876 hypre_ParCSRCommPkgNumRecvs(comm_pkg_S) = num_recvs;
877 send_procs = hypre_ParCSRCommPkgSendProcs(comm_pkg);
878 send_map_starts = hypre_ParCSRCommPkgSendMapStarts(comm_pkg);
879 send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
880 recv_procs = hypre_ParCSRCommPkgRecvProcs(comm_pkg);
881 recv_vec_starts = hypre_ParCSRCommPkgRecvVecStarts(comm_pkg);
882 send_procs_S = NULL;
883 send_map_elmts_S = NULL;
884 if (num_sends)
885 {
886 send_procs_S = hypre_CTAlloc(int,num_sends);
887 send_map_elmts_S = hypre_CTAlloc(int,
888 num_functions*send_map_starts[num_sends]);
889 }
890 send_map_starts_S = hypre_CTAlloc(int,num_sends+1);
891 recv_vec_starts_S = hypre_CTAlloc(int,num_recvs+1);
892 recv_procs_S = NULL;
893 if (num_recvs) recv_procs_S = hypre_CTAlloc(int,num_recvs);
894 send_map_starts_S[0] = 0;
895 for (i=0; i < num_sends; i++)
896 {
897 send_procs_S[i] = send_procs[i];
898 send_map_starts_S[i+1] = num_functions*send_map_starts[i+1];
899 }
900 recv_vec_starts_S[0] = 0;
901 for (i=0; i < num_recvs; i++)
902 {
903 recv_procs_S[i] = recv_procs[i];
904 recv_vec_starts_S[i+1] = num_functions*recv_vec_starts[i+1];
905 }
906
907 cnt = 0;
908 for (i=0; i < send_map_starts[num_sends]; i++)
909 {
910 k1 = num_functions*send_map_elmts[i];
911 for (j=0; j < num_functions; j++)
912 {
913 send_map_elmts_S[cnt++] = k1+j;
914 }
915 }
916 hypre_ParCSRCommPkgSendProcs(comm_pkg_S) = send_procs_S;
917 hypre_ParCSRCommPkgSendMapStarts(comm_pkg_S) = send_map_starts_S;
918 hypre_ParCSRCommPkgSendMapElmts(comm_pkg_S) = send_map_elmts_S;
919 hypre_ParCSRCommPkgRecvProcs(comm_pkg_S) = recv_procs_S;
920 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_S) = recv_vec_starts_S;
921 hypre_ParCSRMatrixCommPkg(S) = comm_pkg_S;
922 }
923
924 if (num_cols_offd_S)
925 {
926 S_offd_j = hypre_CTAlloc(int, S_num_nonzeros_offd);
927 hypre_CSRMatrixJ(S_offd) = S_offd_j;
928
929 col_map_offd_S = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd_S);
930
931 cnt = 0;
932 for (i=0; i < num_cols_offd_SN; i++)
933 {
934 k1 = col_map_offd_SN[i]*num_functions;
935 for (j=0; j < num_functions; j++)
936 col_map_offd_S[cnt++] = k1+j;
937 }
938 hypre_ParCSRMatrixColMapOffd(S) = col_map_offd_S;
939 }
940
941 if (col_offd_SN_to_AN)
942 {
943 col_offd_S_to_A = hypre_CTAlloc(int, num_cols_offd_S);
944
945 cnt = 0;
946 for (i=0; i < num_cols_offd_SN; i++)
947 {
948 k1 = col_offd_SN_to_AN[i]*num_functions;
949 for (j=0; j < num_functions; j++)
950 col_offd_S_to_A[cnt++] = k1+j;
951 }
952 *col_offd_S_to_A_ptr = col_offd_S_to_A;
953 }
954
955
956
957 cnt = 0;
958 row = 0;
959 for (i=0; i < num_nodes; i++)
960 {
961 row++;
962 start = cnt;
963 for (j=SN_diag_i[i]; j < SN_diag_i[i+1]; j++)
964 {
965 jj = SN_diag_j[j];
966 if (data) S_diag_data[cnt] = SN_diag_data[j];
967 S_diag_j[cnt++] = jj*num_functions;
968 }
969 end = cnt;
970 S_diag_i[row] = cnt;
971 for (k1=1; k1 < num_functions; k1++)
972 {
973 row++;
974 for (k=start; k < end; k++)
975 {
976 if (data) S_diag_data[cnt] = S_diag_data[k];
977 S_diag_j[cnt++] = S_diag_j[k]+k1;
978 }
979 S_diag_i[row] = cnt;
980 }
981 }
982
983 cnt = 0;
984 row = 0;
985 for (i=0; i < num_nodes; i++)
986 {
987 row++;
988 start = cnt;
989 for (j=SN_offd_i[i]; j < SN_offd_i[i+1]; j++)
990 {
991 jj = SN_offd_j[j];
992 if (data) S_offd_data[cnt] = SN_offd_data[j];
993 S_offd_j[cnt++] = jj*num_functions;
994 }
995 end = cnt;
996 S_offd_i[row] = cnt;
997 for (k1=1; k1 < num_functions; k1++)
998 {
999 row++;
1000 for (k=start; k < end; k++)
1001 {
1002 if (data) S_offd_data[cnt] = S_offd_data[k];
1003 S_offd_j[cnt++] = S_offd_j[k]+k1;
1004 }
1005 S_offd_i[row] = cnt;
1006 }
1007 }
1008
1009 *S_ptr = S;
1010
1011 return (ierr);
1012}
1013
1014
1015/* This function just finds the scalaer CF_marker and dof_func */
1016
1017int
1018hypre_BoomerAMGCreateScalarCF(int *CFN_marker,
1019 int num_functions,
1020 int num_nodes,
1021 int **dof_func_ptr,
1022 int **CF_marker_ptr)
1023
1024{
1025 int *CF_marker;
1026 int *dof_func;
1027 int num_variables;
1028 int num_coarse_nodes;
1029 int i,j,k,cnt;
1030 int ierr = 0;
1031
1032
1033 num_variables = num_functions*num_nodes;
1034 CF_marker = hypre_CTAlloc(int, num_variables);
1035
1036 cnt = 0;
1037 num_coarse_nodes = 0;
1038 for (i=0; i < num_nodes; i++)
1039 {
1040 if (CFN_marker[i] == 1) num_coarse_nodes++;
1041 for (j=0; j < num_functions; j++)
1042 CF_marker[cnt++] = CFN_marker[i];
1043 }
1044
1045
1046 dof_func = hypre_CTAlloc(int,num_coarse_nodes*num_functions);
1047 cnt = 0;
1048 for (i=0; i < num_nodes; i++)
1049 {
1050 if (CFN_marker[i] == 1)
1051 {
1052 for (k=0; k < num_functions; k++)
1053 dof_func[cnt++] = k;
1054 }
1055 }
1056
1057
1058 *dof_func_ptr = dof_func;
1059 *CF_marker_ptr = CF_marker;
1060
1061
1062 return (ierr);
1063}
Note: See TracBrowser for help on using the repository browser.