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_IJ_mv.h"
9 #include "_hypre_sstruct_ls.h"
10
11 #include "nd1_amge_interpolation.h"
12
13 /*
14 Assume that we are given a fine and coarse topology and the
15 coarse degrees of freedom (DOFs) have been chosen. Assume also,
16 that the global interpolation matrix dof_DOF has a prescribed
17 nonzero pattern. Then, the fine degrees of freedom can be split
18 into 4 groups (here "i" stands for "interior"):
19
20 NODEidof - dofs which are interpolated only from the DOF
21 in one coarse vertex
22 EDGEidof - dofs which are interpolated only from the DOFs
23 in one coarse edge
24 FACEidof - dofs which are interpolated only from the DOFs
25 in one coarse face
26 ELEMidof - dofs which are interpolated only from the DOFs
27 in one coarse element
28
29 The interpolation operator dof_DOF can be build in 4 steps, by
30 consequently filling-in the rows corresponding to the above groups.
31 The code below uses harmonic extension to extend the interpolation
32 from one group to the next.
33 */
hypre_ND1AMGeInterpolation(hypre_ParCSRMatrix * Aee,hypre_ParCSRMatrix * ELEM_idof,hypre_ParCSRMatrix * FACE_idof,hypre_ParCSRMatrix * EDGE_idof,hypre_ParCSRMatrix * ELEM_FACE,hypre_ParCSRMatrix * ELEM_EDGE,HYPRE_Int num_OffProcRows,hypre_MaxwellOffProcRow ** OffProcRows,hypre_IJMatrix * IJ_dof_DOF)34 HYPRE_Int hypre_ND1AMGeInterpolation (hypre_ParCSRMatrix * Aee,
35 hypre_ParCSRMatrix * ELEM_idof,
36 hypre_ParCSRMatrix * FACE_idof,
37 hypre_ParCSRMatrix * EDGE_idof,
38 hypre_ParCSRMatrix * ELEM_FACE,
39 hypre_ParCSRMatrix * ELEM_EDGE,
40 HYPRE_Int num_OffProcRows,
41 hypre_MaxwellOffProcRow ** OffProcRows,
42 hypre_IJMatrix * IJ_dof_DOF)
43 {
44 HYPRE_Int ierr = 0;
45
46 HYPRE_Int i, j;
47 HYPRE_BigInt big_k;
48 HYPRE_BigInt *offproc_rnums;
49 HYPRE_Int *swap;
50
51 hypre_ParCSRMatrix * dof_DOF = (hypre_ParCSRMatrix *)hypre_IJMatrixObject(IJ_dof_DOF);
52 hypre_ParCSRMatrix * ELEM_DOF = ELEM_EDGE;
53 hypre_ParCSRMatrix * ELEM_FACEidof;
54 hypre_ParCSRMatrix * ELEM_EDGEidof;
55 hypre_CSRMatrix *A, *P;
56 HYPRE_Int numELEM = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(ELEM_EDGE));
57
58 HYPRE_Int getrow_ierr;
59 HYPRE_Int three_dimensional_problem;
60
61 MPI_Comm comm= hypre_ParCSRMatrixComm(Aee);
62 HYPRE_Int myproc;
63
64 hypre_MPI_Comm_rank(comm, &myproc);
65
66 #if 0
67 hypre_IJMatrix * ij_dof_DOF = hypre_CTAlloc(hypre_IJMatrix, 1, HYPRE_MEMORY_HOST);
68 /* Convert dof_DOF to IJ matrix, so we can use AddToValues */
69 hypre_IJMatrixComm(ij_dof_DOF) = hypre_ParCSRMatrixComm(dof_DOF);
70 hypre_IJMatrixRowPartitioning(ij_dof_DOF) =
71 hypre_ParCSRMatrixRowStarts(dof_DOF);
72 hypre_IJMatrixColPartitioning(ij_dof_DOF) =
73 hypre_ParCSRMatrixColStarts(dof_DOF);
74 hypre_IJMatrixObject(ij_dof_DOF) = dof_DOF;
75 hypre_IJMatrixAssembleFlag(ij_dof_DOF) = 1;
76 #endif
77
78 /* sort the offproc rows to get quicker comparison for later */
79 if (num_OffProcRows)
80 {
81 offproc_rnums= hypre_TAlloc(HYPRE_BigInt, num_OffProcRows, HYPRE_MEMORY_HOST);
82 swap = hypre_TAlloc(HYPRE_Int, num_OffProcRows, HYPRE_MEMORY_HOST);
83 for (i= 0; i< num_OffProcRows; i++)
84 {
85 offproc_rnums[i]=(OffProcRows[i] -> row);
86 swap[i] = i;
87 }
88 }
89
90 if (num_OffProcRows > 1)
91 {
92 hypre_BigQsortbi(offproc_rnums, swap, 0, num_OffProcRows-1);
93 }
94
95 if (FACE_idof == EDGE_idof)
96 three_dimensional_problem = 0;
97 else
98 three_dimensional_problem = 1;
99
100 /* ELEM_FACEidof = ELEM_FACE x FACE_idof */
101 if (three_dimensional_problem)
102 ELEM_FACEidof = hypre_ParMatmul(ELEM_FACE, FACE_idof);
103
104 /* ELEM_EDGEidof = ELEM_EDGE x EDGE_idof */
105 ELEM_EDGEidof = hypre_ParMatmul(ELEM_EDGE, EDGE_idof);
106
107 /* Loop over local coarse elements */
108 big_k = hypre_ParCSRMatrixFirstRowIndex(ELEM_EDGE);
109 for (i = 0; i < numELEM; i++, big_k++)
110 {
111 HYPRE_Int size1, size2;
112 HYPRE_BigInt *col_ind0, *col_ind1, *col_ind2;
113
114 HYPRE_BigInt *DOF0, *DOF;
115 HYPRE_Int num_DOF;
116 HYPRE_Int num_idof;
117 HYPRE_BigInt *idof0, *idof, *bdof;
118 HYPRE_Int num_bdof;
119
120 HYPRE_Real *boolean_data;
121
122 /* Determine the coarse DOFs */
123 hypre_ParCSRMatrixGetRow (ELEM_DOF, big_k, &num_DOF, &DOF0, &boolean_data);
124 DOF= hypre_TAlloc(HYPRE_BigInt, num_DOF, HYPRE_MEMORY_HOST);
125 for (j= 0; j< num_DOF; j++)
126 {
127 DOF[j]= DOF0[j];
128 }
129 hypre_ParCSRMatrixRestoreRow (ELEM_DOF, big_k, &num_DOF, &DOF0, &boolean_data);
130
131 hypre_BigQsort0(DOF,0,num_DOF-1);
132
133 /* Find the fine dofs interior for the current coarse element */
134 hypre_ParCSRMatrixGetRow (ELEM_idof, big_k, &num_idof, &idof0, &boolean_data);
135 idof= hypre_TAlloc(HYPRE_BigInt, num_idof, HYPRE_MEMORY_HOST);
136 for (j= 0; j< num_idof; j++)
137 {
138 idof[j]= idof0[j];
139 }
140 hypre_ParCSRMatrixRestoreRow (ELEM_idof, big_k, &num_idof, &idof0, &boolean_data);
141
142 /* Sort the interior dofs according to their global number */
143 hypre_BigQsort0(idof,0,num_idof-1);
144
145 /* Find the fine dofs on the boundary of the current coarse element */
146 if (three_dimensional_problem)
147 {
148 hypre_ParCSRMatrixGetRow (ELEM_FACEidof, big_k, &size1, &col_ind0, &boolean_data);
149 col_ind1= hypre_TAlloc(HYPRE_BigInt, size1, HYPRE_MEMORY_HOST);
150 for (j= 0; j< size1; j++)
151 {
152 col_ind1[j]= col_ind0[j];
153 }
154 hypre_ParCSRMatrixRestoreRow (ELEM_FACEidof, big_k, &size1, &col_ind0, &boolean_data);
155 }
156 else
157 size1 = 0;
158
159 hypre_ParCSRMatrixGetRow (ELEM_EDGEidof, big_k, &size2, &col_ind0, &boolean_data);
160 col_ind2= hypre_TAlloc(HYPRE_BigInt, size2, HYPRE_MEMORY_HOST);
161 for (j= 0; j< size2; j++)
162 {
163 col_ind2[j]= col_ind0[j];
164 }
165 hypre_ParCSRMatrixRestoreRow (ELEM_EDGEidof, big_k, &size2, &col_ind0, &boolean_data);
166
167 /* Merge and sort the boundary dofs according to their global number */
168 num_bdof = size1 + size2;
169 bdof = hypre_CTAlloc(HYPRE_BigInt, num_bdof, HYPRE_MEMORY_HOST);
170 if (three_dimensional_problem)
171 hypre_TMemcpy(bdof, col_ind1, HYPRE_BigInt, size1, HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
172 hypre_TMemcpy(bdof+size1, col_ind2, HYPRE_BigInt, size2, HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
173
174 hypre_BigQsort0(bdof,0,num_bdof-1);
175
176 /* A = extract_rows(Aee, idof) */
177 A = hypre_CSRMatrixCreate (num_idof, num_idof + num_bdof,
178 num_idof * (num_idof + num_bdof));
179 hypre_CSRMatrixBigInitialize(A);
180 {
181 HYPRE_Int *I = hypre_CSRMatrixI(A);
182 HYPRE_BigInt *J = hypre_CSRMatrixBigJ(A);
183 HYPRE_Real *data = hypre_CSRMatrixData(A);
184
185 HYPRE_BigInt *tmp_J;
186 HYPRE_Real *tmp_data;
187
188 I[0] = 0;
189 for (j = 0; j < num_idof; j++)
190 {
191 getrow_ierr= hypre_ParCSRMatrixGetRow (Aee, idof[j], &I[j+1], &tmp_J, &tmp_data);
192 if (getrow_ierr <0)
193 hypre_printf("getrow Aee off proc[%d] = \n",myproc);
194 hypre_TMemcpy(J, tmp_J, HYPRE_BigInt, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
195 hypre_TMemcpy(data, tmp_data, HYPRE_Real, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
196 J+= I[j+1];
197 data+= I[j+1];
198 hypre_ParCSRMatrixRestoreRow (Aee, idof[j], &I[j+1], &tmp_J, &tmp_data);
199 I[j+1] += I[j];
200 }
201 }
202
203 /* P = extract_rows(dof_DOF, idof+bdof) */
204 P = hypre_CSRMatrixCreate (num_idof + num_bdof, num_DOF,
205 (num_idof + num_bdof) * num_DOF);
206 hypre_CSRMatrixBigInitialize(P);
207 {
208 HYPRE_Int *I = hypre_CSRMatrixI(P);
209 HYPRE_BigInt *J = hypre_CSRMatrixBigJ(P);
210 HYPRE_Real *data = hypre_CSRMatrixData(P);
211 HYPRE_Int m;
212
213 HYPRE_BigInt *tmp_J;
214 HYPRE_Real *tmp_data;
215
216 I[0] = 0;
217 for (j = 0; j < num_idof; j++)
218 {
219 getrow_ierr= hypre_ParCSRMatrixGetRow (dof_DOF, idof[j], &I[j+1], &tmp_J, &tmp_data);
220 if (getrow_ierr >= 0)
221 {
222 hypre_TMemcpy(J, tmp_J, HYPRE_BigInt, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
223 hypre_TMemcpy(data, tmp_data, HYPRE_Real, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
224 J+= I[j+1];
225 data+= I[j+1];
226 hypre_ParCSRMatrixRestoreRow (dof_DOF, idof[j], &I[j+1], &tmp_J, &tmp_data);
227 I[j+1] += I[j];
228 }
229 else /* row offproc */
230 {
231 hypre_ParCSRMatrixRestoreRow (dof_DOF, idof[j], &I[j+1], &tmp_J, &tmp_data);
232 /* search for OffProcRows */
233 m= 0;
234 while (m < num_OffProcRows)
235 {
236 if (offproc_rnums[m] == idof[j])
237 {
238 break;
239 }
240 else
241 {
242 m++;
243 }
244 }
245 I[j+1]= (OffProcRows[swap[m]] -> ncols);
246 tmp_J = (OffProcRows[swap[m]] -> cols);
247 tmp_data= (OffProcRows[swap[m]] -> data);
248 hypre_TMemcpy(J, tmp_J, HYPRE_BigInt, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
249 hypre_TMemcpy(data, tmp_data, HYPRE_Real, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
250 J+= I[j+1];
251 data+= I[j+1];
252 I[j+1] += I[j];
253 }
254
255 }
256 for ( ; j < num_idof + num_bdof; j++)
257 {
258 getrow_ierr= hypre_ParCSRMatrixGetRow (dof_DOF, bdof[j-num_idof], &I[j+1], &tmp_J, &tmp_data);
259 if (getrow_ierr >= 0)
260 {
261 hypre_TMemcpy(J, tmp_J, HYPRE_BigInt, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
262 hypre_TMemcpy(data, tmp_data, HYPRE_Real, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
263 J+= I[j+1];
264 data+= I[j+1];
265 hypre_ParCSRMatrixRestoreRow (dof_DOF, bdof[j-num_idof], &I[j+1], &tmp_J, &tmp_data);
266 I[j+1] += I[j];
267 }
268 else /* row offproc */
269 {
270 hypre_ParCSRMatrixRestoreRow (dof_DOF, bdof[j-num_idof], &I[j+1], &tmp_J, &tmp_data);
271 /* search for OffProcRows */
272 m= 0;
273 while (m < num_OffProcRows)
274 {
275 if (offproc_rnums[m] == bdof[j-num_idof])
276 {
277 break;
278 }
279 else
280 {
281 m++;
282 }
283 }
284 if (m>= num_OffProcRows)hypre_printf("here the mistake\n");
285 I[j+1]= (OffProcRows[swap[m]] -> ncols);
286 tmp_J = (OffProcRows[swap[m]] -> cols);
287 tmp_data= (OffProcRows[swap[m]] -> data);
288 hypre_TMemcpy(J, tmp_J, HYPRE_BigInt, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
289 hypre_TMemcpy(data, tmp_data, HYPRE_Real, I[j+1], HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
290 J+= I[j+1];
291 data+= I[j+1];
292 I[j+1] += I[j];
293 }
294 }
295 }
296
297 /* Pi = Aii^{-1} Aib Pb */
298 hypre_HarmonicExtension (A, P, num_DOF, DOF,
299 num_idof, idof, num_bdof, bdof);
300
301 /* Insert Pi in dof_DOF */
302 {
303 HYPRE_Int * ncols = hypre_CTAlloc(HYPRE_Int, num_idof, HYPRE_MEMORY_HOST);
304 HYPRE_Int * idof_indexes = hypre_CTAlloc(HYPRE_Int, num_idof, HYPRE_MEMORY_HOST);
305
306 for (j = 0; j < num_idof; j++)
307 {
308 ncols[j] = num_DOF;
309 idof_indexes[j] = j*num_DOF;
310 }
311
312 hypre_IJMatrixAddToValuesParCSR (IJ_dof_DOF,
313 num_idof, ncols, idof, idof_indexes,
314 hypre_CSRMatrixBigJ(P),
315 hypre_CSRMatrixData(P));
316
317 hypre_TFree(ncols, HYPRE_MEMORY_HOST);
318 hypre_TFree(idof_indexes, HYPRE_MEMORY_HOST);
319 }
320
321 hypre_TFree(DOF, HYPRE_MEMORY_HOST);
322 hypre_TFree(idof, HYPRE_MEMORY_HOST);
323 if (three_dimensional_problem)
324 {
325 hypre_TFree(col_ind1, HYPRE_MEMORY_HOST);
326 }
327 hypre_TFree(col_ind2, HYPRE_MEMORY_HOST);
328 hypre_TFree(bdof, HYPRE_MEMORY_HOST);
329
330 hypre_CSRMatrixDestroy(A);
331 hypre_CSRMatrixDestroy(P);
332 }
333
334 #if 0
335 hypre_TFree(ij_dof_DOF, HYPRE_MEMORY_HOST);
336 #endif
337
338 if (three_dimensional_problem)
339 hypre_ParCSRMatrixDestroy(ELEM_FACEidof);
340 hypre_ParCSRMatrixDestroy(ELEM_EDGEidof);
341
342 if (num_OffProcRows)
343 {
344 hypre_TFree(offproc_rnums, HYPRE_MEMORY_HOST);
345 hypre_TFree(swap, HYPRE_MEMORY_HOST);
346 }
347
348 return ierr;
349 }
350
351
352
353
hypre_HarmonicExtension(hypre_CSRMatrix * A,hypre_CSRMatrix * P,HYPRE_Int num_DOF,HYPRE_BigInt * DOF,HYPRE_Int num_idof,HYPRE_BigInt * idof,HYPRE_Int num_bdof,HYPRE_BigInt * bdof)354 HYPRE_Int hypre_HarmonicExtension (hypre_CSRMatrix *A,
355 hypre_CSRMatrix *P,
356 HYPRE_Int num_DOF, HYPRE_BigInt *DOF,
357 HYPRE_Int num_idof, HYPRE_BigInt *idof,
358 HYPRE_Int num_bdof, HYPRE_BigInt *bdof)
359 {
360 HYPRE_Int ierr = 0;
361
362 HYPRE_Int i, j, k, l, m;
363 HYPRE_Real factor;
364
365 HYPRE_Int *IA = hypre_CSRMatrixI(A);
366 HYPRE_BigInt *JA = hypre_CSRMatrixBigJ(A);
367 HYPRE_Real *dataA = hypre_CSRMatrixData(A);
368
369 HYPRE_Int *IP = hypre_CSRMatrixI(P);
370 HYPRE_BigInt *JP = hypre_CSRMatrixBigJ(P);
371 HYPRE_Real *dataP = hypre_CSRMatrixData(P);
372
373 HYPRE_Real * Aii = hypre_CTAlloc(HYPRE_Real, num_idof*num_idof, HYPRE_MEMORY_HOST);
374 HYPRE_Real * Pi = hypre_CTAlloc(HYPRE_Real, num_idof*num_DOF, HYPRE_MEMORY_HOST);
375
376 /* Loop over the rows of A */
377 for (i = 0; i < num_idof; i++)
378 for (j = IA[i]; j < IA[i+1]; j++)
379 {
380 /* Global to local*/
381 k = hypre_BigBinarySearch(idof,JA[j], num_idof);
382 /* If a column is a bdof, compute its participation in Pi = Aib x Pb */
383 if (k == -1)
384 {
385 k = hypre_BigBinarySearch(bdof,JA[j], num_bdof);
386 if (k > -1)
387 {
388 for (l = IP[k+num_idof]; l < IP[k+num_idof+1]; l++)
389 {
390 m = hypre_BigBinarySearch(DOF,JP[l], num_DOF);
391 if (m > -1)
392 {
393 m+=i*num_DOF;
394 /* Pi[i*num_DOF+m] += dataA[j] * dataP[l];*/
395 Pi[m] += dataA[j] * dataP[l];
396 }
397 }
398 }
399 }
400 /* If a column is an idof, put it in Aii */
401 else
402 Aii[i*num_idof+k] = dataA[j];
403 }
404
405 /* Perform Gaussian elimination in [Aii, Pi] */
406 for (j = 0; j < num_idof-1; j++)
407 if (Aii[j*num_idof+j] != 0.0)
408 for (i = j+1; i < num_idof; i++)
409 if (Aii[i*num_idof+j] != 0.0)
410 {
411 factor = Aii[i*num_idof+j]/Aii[j*num_idof+j];
412 for (m = j+1; m < num_idof; m++)
413 Aii[i*num_idof+m] -= factor * Aii[j*num_idof+m];
414 for (m = 0; m < num_DOF; m++)
415 Pi[i*num_DOF+m] -= factor * Pi[j*num_DOF+m];
416 }
417
418 /* Back Substitution */
419 for (i = num_idof-1; i >= 0; i--)
420 {
421 for (j = i+1; j < num_idof; j++)
422 if (Aii[i*num_idof+j] != 0.0)
423 for (m = 0; m < num_DOF; m++)
424 Pi[i*num_DOF+m] -= Aii[i*num_idof+j] * Pi[j*num_DOF+m];
425
426 for (m = 0; m < num_DOF; m++)
427 Pi[i*num_DOF+m] /= Aii[i*num_idof+i];
428 }
429
430 /* Put -Pi back in P. We assume that each idof depends on _all_ DOFs */
431 for (i = 0; i < num_idof; i++, JP += num_DOF, dataP += num_DOF)
432 for (j = 0; j < num_DOF; j++)
433 {
434 JP[j] = DOF[j];
435 dataP[j] = -Pi[i*num_DOF+j];
436 }
437
438 hypre_TFree(Aii, HYPRE_MEMORY_HOST);
439 hypre_TFree(Pi, HYPRE_MEMORY_HOST);
440
441 return ierr;
442 }
443