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 /******************************************************************************
9 *
10 * DiagScale - Diagonal scaling.
11 *
12 *****************************************************************************/
13
14 #include <stdlib.h>
15 #include "math.h"
16 #include "Common.h"
17 #include "Matrix.h"
18 #include "RowPatt.h"
19 #include "DiagScale.h"
20 #include "OrderStat.h"
21 #include "Mem.h"
22
23 HYPRE_Int FindNumReplies(MPI_Comm comm, HYPRE_Int *replies_list);
24
25 #define DIAG_VALS_TAG 225
26 #define DIAG_INDS_TAG 226
27
28 /*--------------------------------------------------------------------------
29 * ExchangeDiagEntries - Given a list of indices of diagonal entries required
30 * by this processor, "reqind" of length "reqlen", return a list of
31 * corresponding diagonal entries, "diags". Used internally only by
32 * DiagScaleCreate.
33 *
34 * comm - MPI communicator (input)
35 * mat - matrix used to map row and column numbers to processors (input)
36 * reqlen - length of request list (input)
37 * reqind - list of indices (input)
38 * diags - corresponding list of diagonal entries (output)
39 * num_requests - number of requests (output)
40 * requests - request handles, used to check that all responses are back
41 * (output)
42 * replies_list - array that indicates who we sent message to (output)
43 *--------------------------------------------------------------------------*/
44
ExchangeDiagEntries(MPI_Comm comm,Matrix * mat,HYPRE_Int reqlen,HYPRE_Int * reqind,HYPRE_Real * diags,HYPRE_Int * num_requests,hypre_MPI_Request * requests,HYPRE_Int * replies_list)45 static void ExchangeDiagEntries(MPI_Comm comm, Matrix *mat, HYPRE_Int reqlen,
46 HYPRE_Int *reqind, HYPRE_Real *diags, HYPRE_Int *num_requests, hypre_MPI_Request *requests,
47 HYPRE_Int *replies_list)
48 {
49 hypre_MPI_Request request;
50 HYPRE_Int i, j, this_pe;
51
52 hypre_shell_sort(reqlen, reqind);
53
54 *num_requests = 0;
55
56 for (i=0; i<reqlen; i=j) /* j is set below */
57 {
58 /* The processor that owns the row with index reqind[i] */
59 this_pe = MatrixRowPe(mat, reqind[i]);
60
61 /* Figure out other rows we need from this_pe */
62 for (j=i+1; j<reqlen; j++)
63 {
64 /* if row is on different pe */
65 if (reqind[j] < mat->beg_rows[this_pe] ||
66 reqind[j] > mat->end_rows[this_pe])
67 break;
68 }
69
70 /* Post receive for diagonal values */
71 hypre_MPI_Irecv(&diags[i], j-i, hypre_MPI_REAL, this_pe, DIAG_VALS_TAG,
72 comm, &requests[*num_requests]);
73
74 /* Request rows in reqind[i..j-1] */
75 hypre_MPI_Isend(&reqind[i], j-i, HYPRE_MPI_INT, this_pe, DIAG_INDS_TAG,
76 comm, &request);
77 hypre_MPI_Request_free(&request);
78 (*num_requests)++;
79
80 if (replies_list != NULL)
81 replies_list[this_pe] = 1;
82 }
83 }
84
85 /*--------------------------------------------------------------------------
86 * ExchangeDiagEntriesServer - Receive requests for diagonal entries and
87 * send replies. Used internally only by DiagScaleCreate.
88 *
89 * comm - MPI communicator (input)
90 * mat - matrix used to map row and column numbers to processors (input)
91 * local_diags - local diagonal entries (input)
92 * num_requests - number of requests to be received (input)
93 *--------------------------------------------------------------------------*/
94
ExchangeDiagEntriesServer(MPI_Comm comm,Matrix * mat,HYPRE_Real * local_diags,HYPRE_Int num_requests,Mem * mem,hypre_MPI_Request * requests)95 static void ExchangeDiagEntriesServer(MPI_Comm comm, Matrix *mat,
96 HYPRE_Real *local_diags, HYPRE_Int num_requests, Mem *mem, hypre_MPI_Request *requests)
97 {
98 hypre_MPI_Status status;
99 HYPRE_Int *recvbuf;
100 HYPRE_Real *sendbuf;
101 HYPRE_Int i, j, source, count;
102
103 /* recvbuf contains requested indices */
104 /* sendbuf contains corresponding diagonal entries */
105
106 for (i=0; i<num_requests; i++)
107 {
108 hypre_MPI_Probe(hypre_MPI_ANY_SOURCE, DIAG_INDS_TAG, comm, &status);
109 source = status.hypre_MPI_SOURCE;
110 hypre_MPI_Get_count(&status, HYPRE_MPI_INT, &count);
111
112 recvbuf = (HYPRE_Int *) MemAlloc(mem, count*sizeof(HYPRE_Int));
113 sendbuf = (HYPRE_Real *) MemAlloc(mem, count*sizeof(HYPRE_Real));
114
115 /*hypre_MPI_Recv(recvbuf, count, HYPRE_MPI_INT, hypre_MPI_ANY_SOURCE, */
116 hypre_MPI_Recv(recvbuf, count, HYPRE_MPI_INT, source,
117 DIAG_INDS_TAG, comm, &status);
118 source = status.hypre_MPI_SOURCE;
119
120 /* Construct reply message of diagonal entries in sendbuf */
121 for (j=0; j<count; j++)
122 sendbuf[j] = local_diags[recvbuf[j] - mat->beg_row];
123
124 /* Use ready-mode send, since receives already posted */
125 hypre_MPI_Irsend(sendbuf, count, hypre_MPI_REAL, source,
126 DIAG_VALS_TAG, comm, &requests[i]);
127 }
128 }
129
130 /*--------------------------------------------------------------------------
131 * DiagScaleCreate - Return (a pointer to) a diagonal scaling object.
132 * Scale using the diagonal of A. Use the list of external indices
133 * from the numbering object "numb".
134 *--------------------------------------------------------------------------*/
135
DiagScaleCreate(Matrix * A,Numbering * numb)136 DiagScale *DiagScaleCreate(Matrix *A, Numbering *numb)
137 {
138 hypre_MPI_Request *requests;
139 hypre_MPI_Status *statuses;
140 HYPRE_Int npes, row, j, num_requests, num_replies, *replies_list;
141 HYPRE_Int len, *ind;
142 HYPRE_Real *val, *temp;
143
144 Mem *mem;
145 hypre_MPI_Request *requests2;
146
147 DiagScale *p = hypre_TAlloc(DiagScale, 1, HYPRE_MEMORY_HOST);
148
149 /* Storage for local diagonal entries */
150 p->local_diags = (HYPRE_Real *)
151 hypre_TAlloc(HYPRE_Real, (A->end_row - A->beg_row + 1) , HYPRE_MEMORY_HOST);
152
153 /* Extract the local diagonal entries */
154 for (row=0; row<=A->end_row - A->beg_row; row++)
155 {
156 MatrixGetRow(A, row, &len, &ind, &val);
157
158 p->local_diags[row] = 1.0; /* in case no diag entry */
159
160 for (j=0; j<len; j++)
161 {
162 if (ind[j] == row)
163 {
164 if (val[j] != 0.0)
165 p->local_diags[row] = 1.0 / sqrt(ABS(val[j]));
166 break;
167 }
168 }
169 }
170
171 /* Get the list of diagonal indices that we need.
172 This is simply the external indices */
173 /* ExchangeDiagEntries will sort the list - so give it a copy */
174 len = numb->num_ind - numb->num_loc;
175 ind = NULL;
176 p->ext_diags = NULL;
177 if (len)
178 {
179 ind = hypre_TAlloc(HYPRE_Int, len , HYPRE_MEMORY_HOST);
180 hypre_TMemcpy(ind, &numb->local_to_global[numb->num_loc], HYPRE_Int, len, HYPRE_MEMORY_HOST, HYPRE_MEMORY_HOST);
181
182 /* buffer for receiving diagonal values from other processors */
183 p->ext_diags = hypre_TAlloc(HYPRE_Real, len , HYPRE_MEMORY_HOST);
184 }
185
186 hypre_MPI_Comm_size(A->comm, &npes);
187 requests = hypre_TAlloc(hypre_MPI_Request, npes , HYPRE_MEMORY_HOST);
188 statuses = hypre_TAlloc(hypre_MPI_Status, npes , HYPRE_MEMORY_HOST);
189 replies_list = hypre_CTAlloc(HYPRE_Int, npes, HYPRE_MEMORY_HOST);
190
191 ExchangeDiagEntries(A->comm, A, len, ind, p->ext_diags, &num_requests,
192 requests, replies_list);
193
194 num_replies = FindNumReplies(A->comm, replies_list);
195 hypre_TFree(replies_list,HYPRE_MEMORY_HOST);
196
197 mem = MemCreate();
198 requests2 = NULL;
199 if (num_replies)
200 requests2 = hypre_TAlloc(hypre_MPI_Request, num_replies , HYPRE_MEMORY_HOST);
201
202 ExchangeDiagEntriesServer(A->comm, A, p->local_diags, num_replies,
203 mem, requests2);
204
205 /* Wait for all replies */
206 hypre_MPI_Waitall(num_requests, requests, statuses);
207 hypre_TFree(requests,HYPRE_MEMORY_HOST);
208
209 p->offset = A->end_row - A->beg_row + 1;
210
211 /* ind contains global indices corresponding to order that entries
212 are stored in ext_diags. Reorder ext_diags in original ordering */
213 NumberingGlobalToLocal(numb, len, ind, ind);
214 temp = NULL;
215 if (len)
216 temp = hypre_TAlloc(HYPRE_Real, len , HYPRE_MEMORY_HOST);
217 for (j=0; j<len; j++)
218 temp[ind[j]-p->offset] = p->ext_diags[j];
219
220 hypre_TFree(ind,HYPRE_MEMORY_HOST);
221 hypre_TFree(p->ext_diags,HYPRE_MEMORY_HOST);
222 p->ext_diags = temp;
223
224 /* Wait for all sends */
225 hypre_MPI_Waitall(num_replies, requests2, statuses);
226 hypre_TFree(requests2,HYPRE_MEMORY_HOST);
227 MemDestroy(mem);
228
229 hypre_TFree(statuses,HYPRE_MEMORY_HOST);
230 return p;
231 }
232
233 /*--------------------------------------------------------------------------
234 * DiagScaleDestroy - Destroy a diagonal scale object.
235 *--------------------------------------------------------------------------*/
236
DiagScaleDestroy(DiagScale * p)237 void DiagScaleDestroy(DiagScale *p)
238 {
239 hypre_TFree(p->local_diags,HYPRE_MEMORY_HOST);
240 hypre_TFree(p->ext_diags,HYPRE_MEMORY_HOST);
241
242 hypre_TFree(p,HYPRE_MEMORY_HOST);
243 }
244
245 /*--------------------------------------------------------------------------
246 * DiagScaleGet - Returns scale factor given a row number in local indexing.
247 * The factor is the reciprocal of the square root of the diagonal entry.
248 *--------------------------------------------------------------------------*/
249
DiagScaleGet(DiagScale * p,HYPRE_Int index)250 HYPRE_Real DiagScaleGet(DiagScale *p, HYPRE_Int index)
251 {
252 if (index < p->offset)
253 {
254 return p->local_diags[index];
255 }
256 else
257 {
258 return p->ext_diags[index - p->offset];
259 }
260 }
261