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_parcsr_block_mv.h"
9
10 /*---------------------------------------------------------------------------
11 * hypre_BoomerAMGBlockCreateNodalA
12
13 This is the block version of creating a nodal norm matrix.
14
15 option: determine which type of "norm" (or other measurement) is used.
16
17 1 = frobenius
18 2 = sum of abs. value of all elements
19 3 = largest element (positive or negative)
20 4 = 1-norm
21 5 = inf - norm
22 6 = sum of all elements
23 *--------------------------------------------------------------------------*/
24
25 HYPRE_Int
hypre_BoomerAMGBlockCreateNodalA(hypre_ParCSRBlockMatrix * A,HYPRE_Int option,HYPRE_Int diag_option,hypre_ParCSRMatrix ** AN_ptr)26 hypre_BoomerAMGBlockCreateNodalA(hypre_ParCSRBlockMatrix *A,
27 HYPRE_Int option,
28 HYPRE_Int diag_option,
29 hypre_ParCSRMatrix **AN_ptr)
30 {
31 MPI_Comm comm = hypre_ParCSRBlockMatrixComm(A);
32 hypre_CSRBlockMatrix *A_diag = hypre_ParCSRBlockMatrixDiag(A);
33 HYPRE_Int *A_diag_i = hypre_CSRBlockMatrixI(A_diag);
34 HYPRE_Real *A_diag_data = hypre_CSRBlockMatrixData(A_diag);
35
36 HYPRE_Int block_size = hypre_CSRBlockMatrixBlockSize(A_diag);
37 HYPRE_Int bnnz = block_size*block_size;
38
39 hypre_CSRBlockMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
40 HYPRE_Int *A_offd_i = hypre_CSRBlockMatrixI(A_offd);
41 HYPRE_Real *A_offd_data = hypre_CSRBlockMatrixData(A_offd);
42 HYPRE_Int *A_diag_j = hypre_CSRBlockMatrixJ(A_diag);
43 HYPRE_Int *A_offd_j = hypre_CSRBlockMatrixJ(A_offd);
44
45 HYPRE_BigInt *row_starts = hypre_ParCSRBlockMatrixRowStarts(A);
46 HYPRE_BigInt *col_map_offd = hypre_ParCSRBlockMatrixColMapOffd(A);
47 HYPRE_Int num_nonzeros_diag;
48 HYPRE_Int num_nonzeros_offd = 0;
49 HYPRE_Int num_cols_offd = 0;
50
51 hypre_ParCSRMatrix *AN;
52 hypre_CSRMatrix *AN_diag;
53 HYPRE_Int *AN_diag_i;
54 HYPRE_Int *AN_diag_j=NULL;
55 HYPRE_Real *AN_diag_data = NULL;
56 hypre_CSRMatrix *AN_offd;
57 HYPRE_Int *AN_offd_i;
58 HYPRE_Int *AN_offd_j = NULL;
59 HYPRE_Real *AN_offd_data = NULL;
60 HYPRE_BigInt *col_map_offd_AN = NULL;
61
62 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRBlockMatrixCommPkg(A);
63 HYPRE_Int num_sends;
64 HYPRE_Int num_recvs;
65 HYPRE_Int *send_procs;
66 HYPRE_Int *send_map_starts;
67 HYPRE_Int *send_map_elmts;
68 HYPRE_Int *recv_procs;
69 HYPRE_Int *recv_vec_starts;
70
71 hypre_ParCSRCommPkg *comm_pkg_AN = NULL;
72 HYPRE_Int *send_procs_AN = NULL;
73 HYPRE_Int *send_map_starts_AN = NULL;
74 HYPRE_Int *send_map_elmts_AN = NULL;
75 HYPRE_Int *recv_procs_AN = NULL;
76 HYPRE_Int *recv_vec_starts_AN = NULL;
77
78 HYPRE_Int i;
79
80 HYPRE_Int num_procs;
81 HYPRE_Int cnt;
82 HYPRE_Int norm_type;
83
84 HYPRE_BigInt global_num_nodes;
85 HYPRE_Int num_nodes;
86
87 HYPRE_Int index, k;
88
89 HYPRE_Real tmp;
90 HYPRE_Real sum;
91
92 hypre_MPI_Comm_size(comm,&num_procs);
93
94 if (!comm_pkg)
95 {
96 hypre_BlockMatvecCommPkgCreate(A);
97 comm_pkg = hypre_ParCSRBlockMatrixCommPkg(A);
98 }
99
100 norm_type = hypre_abs(option);
101
102 /* Set up the new matrix AN */
103
104 global_num_nodes = hypre_ParCSRBlockMatrixGlobalNumRows(A);
105 num_nodes = hypre_CSRBlockMatrixNumRows(A_diag);
106
107 /* the diag part */
108
109 num_nonzeros_diag = A_diag_i[num_nodes];
110 AN_diag_i = hypre_CTAlloc(HYPRE_Int, num_nodes+1, HYPRE_MEMORY_HOST);
111
112 for (i = 0; i <= num_nodes; i++)
113 {
114 AN_diag_i[i] = A_diag_i[i];
115 }
116
117 AN_diag_j = hypre_CTAlloc(HYPRE_Int, num_nonzeros_diag, HYPRE_MEMORY_HOST);
118 AN_diag_data = hypre_CTAlloc(HYPRE_Real, num_nonzeros_diag, HYPRE_MEMORY_HOST);
119
120 AN_diag = hypre_CSRMatrixCreate(num_nodes, num_nodes, num_nonzeros_diag);
121 hypre_CSRMatrixI(AN_diag) = AN_diag_i;
122 hypre_CSRMatrixJ(AN_diag) = AN_diag_j;
123 hypre_CSRMatrixData(AN_diag) = AN_diag_data;
124
125 for (i = 0; i < num_nonzeros_diag; i++)
126 {
127 AN_diag_j[i] = A_diag_j[i];
128 hypre_CSRBlockMatrixBlockNorm(norm_type, &A_diag_data[i*bnnz],
129 &tmp, block_size);
130 AN_diag_data[i] = tmp;
131 }
132
133
134 if (diag_option == 1)
135 {
136 /* make the diag entry the negative of the sum of off-diag entries (NEED
137 * to get more below!)*/
138 /* the diagonal is the first element listed in each row - */
139 for (i = 0; i < num_nodes; i++)
140 {
141 index = AN_diag_i[i];
142 sum = 0.0;
143 for (k = AN_diag_i[i]+1; k < AN_diag_i[i+1]; k++)
144 {
145 sum += AN_diag_data[k];
146
147 }
148
149 AN_diag_data[index] = -sum;
150 }
151 }
152 else if (diag_option == 2)
153 {
154 /* make all diagonal entries negative */
155 /* the diagonal is the first element listed in each row - */
156
157 for (i = 0; i < num_nodes; i++)
158 {
159 index = AN_diag_i[i];
160 AN_diag_data[index] = -AN_diag_data[index];
161 }
162 }
163
164 /* copy the commpkg */
165 if (comm_pkg)
166 {
167 comm_pkg_AN = hypre_CTAlloc(hypre_ParCSRCommPkg, 1, HYPRE_MEMORY_HOST);
168 hypre_ParCSRCommPkgComm(comm_pkg_AN) = comm;
169
170 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
171 hypre_ParCSRCommPkgNumSends(comm_pkg_AN) = num_sends;
172
173 num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg);
174 hypre_ParCSRCommPkgNumRecvs(comm_pkg_AN) = num_recvs;
175
176 send_procs = hypre_ParCSRCommPkgSendProcs(comm_pkg);
177 send_map_starts = hypre_ParCSRCommPkgSendMapStarts(comm_pkg);
178 send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
179 if (num_sends)
180 {
181 send_procs_AN = hypre_CTAlloc(HYPRE_Int, num_sends, HYPRE_MEMORY_HOST);
182 send_map_elmts_AN = hypre_CTAlloc(HYPRE_Int, send_map_starts[num_sends], HYPRE_MEMORY_HOST);
183 }
184 send_map_starts_AN = hypre_CTAlloc(HYPRE_Int, num_sends+1, HYPRE_MEMORY_HOST);
185 send_map_starts_AN[0] = 0;
186 for (i = 0; i < num_sends; i++)
187 {
188 send_procs_AN[i] = send_procs[i];
189 send_map_starts_AN[i+1] = send_map_starts[i+1];
190 }
191 cnt = send_map_starts_AN[num_sends];
192 for (i = 0; i < cnt; i++)
193 {
194 send_map_elmts_AN[i] = send_map_elmts[i];
195 }
196 hypre_ParCSRCommPkgSendProcs(comm_pkg_AN) = send_procs_AN;
197 hypre_ParCSRCommPkgSendMapStarts(comm_pkg_AN) = send_map_starts_AN;
198 hypre_ParCSRCommPkgSendMapElmts(comm_pkg_AN) = send_map_elmts_AN;
199
200 recv_procs = hypre_ParCSRCommPkgRecvProcs(comm_pkg);
201 recv_vec_starts = hypre_ParCSRCommPkgRecvVecStarts(comm_pkg);
202 recv_vec_starts_AN = hypre_CTAlloc(HYPRE_Int, num_recvs+1, HYPRE_MEMORY_HOST);
203 if (num_recvs) recv_procs_AN = hypre_CTAlloc(HYPRE_Int, num_recvs, HYPRE_MEMORY_HOST);
204
205 recv_vec_starts_AN[0] = recv_vec_starts[0];
206 for (i = 0; i < num_recvs; i++)
207 {
208 recv_procs_AN[i] = recv_procs[i];
209 recv_vec_starts_AN[i+1] = recv_vec_starts[i+1];
210
211 }
212 hypre_ParCSRCommPkgRecvProcs(comm_pkg_AN) = recv_procs_AN;
213 hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_AN) = recv_vec_starts_AN;
214
215 }
216
217 /* the off-diag part */
218
219 num_cols_offd = hypre_CSRBlockMatrixNumCols(A_offd);
220 col_map_offd_AN = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd, HYPRE_MEMORY_HOST);
221 for (i = 0; i < num_cols_offd; i++)
222 {
223 col_map_offd_AN[i] = col_map_offd[i];
224 }
225
226 num_nonzeros_offd = A_offd_i[num_nodes];
227 AN_offd_i = hypre_CTAlloc(HYPRE_Int, num_nodes+1, HYPRE_MEMORY_HOST);
228 for (i = 0; i <= num_nodes; i++)
229 {
230 AN_offd_i[i] = A_offd_i[i];
231 }
232
233 AN_offd_j = hypre_CTAlloc(HYPRE_Int, num_nonzeros_offd, HYPRE_MEMORY_HOST);
234 AN_offd_data = hypre_CTAlloc(HYPRE_Real, num_nonzeros_offd, HYPRE_MEMORY_HOST);
235
236 for (i = 0; i < num_nonzeros_offd; i++)
237 {
238 AN_offd_j[i] = A_offd_j[i];
239 hypre_CSRBlockMatrixBlockNorm(norm_type, &A_offd_data[i*bnnz],
240 &tmp, block_size);
241 AN_offd_data[i] = tmp;
242 }
243
244 AN_offd = hypre_CSRMatrixCreate(num_nodes, num_cols_offd, num_nonzeros_offd);
245
246 hypre_CSRMatrixI(AN_offd) = AN_offd_i;
247 hypre_CSRMatrixJ(AN_offd) = AN_offd_j;
248 hypre_CSRMatrixData(AN_offd) = AN_offd_data;
249
250 if (diag_option == 1)
251 {
252 /* make the diag entry the negative of the sum of off-diag entries (here
253 we are adding the off_diag contribution)*/
254 /* the diagonal is the first element listed in each row of AN_diag_data - */
255 for (i = 0; i < num_nodes; i++)
256 {
257 sum = 0.0;
258 for (k = AN_offd_i[i]; k < AN_offd_i[i+1]; k++)
259 {
260 sum += AN_offd_data[k];
261 }
262 index = AN_diag_i[i];/* location of diag entry in data */
263 AN_diag_data[index] -= sum; /* subtract from current value */
264 }
265 }
266
267 /* now create AN */
268 AN = hypre_ParCSRMatrixCreate(comm, global_num_nodes, global_num_nodes,
269 row_starts, row_starts, num_cols_offd,
270 num_nonzeros_diag, num_nonzeros_offd);
271
272 /* we already created the diag and offd matrices - so we don't need the ones
273 created above */
274 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(AN));
275 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(AN));
276 hypre_ParCSRMatrixDiag(AN) = AN_diag;
277 hypre_ParCSRMatrixOffd(AN) = AN_offd;
278
279 hypre_CSRMatrixMemoryLocation(AN_diag) = HYPRE_MEMORY_HOST;
280 hypre_CSRMatrixMemoryLocation(AN_offd) = HYPRE_MEMORY_HOST;
281
282 hypre_ParCSRMatrixColMapOffd(AN) = col_map_offd_AN;
283 hypre_ParCSRMatrixCommPkg(AN) = comm_pkg_AN;
284
285 *AN_ptr = AN;
286
287 return hypre_error_flag;
288 }
289