Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
ilu_mpi_bj.c
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43/* to do: re-integrate fix-smalll-pivots */
44
45#include "ilu_dh.h"
46#include "Mem_dh.h"
47#include "Parser_dh.h"
48#include "Euclid_dh.h"
49#include "getRow_dh.h"
50#include "Factor_dh.h"
51#include "SubdomainGraph_dh.h"
52
53
54int symbolic_row_private (int localRow, int beg_row, int end_row,
55 int *list, int *marker, int *tmpFill,
56 int len, int *CVAL, double *AVAL,
57 int *o2n_col, Euclid_dh ctx);
58
59static int numeric_row_private (int localRow, int beg_row, int end_row,
60 int len, int *CVAL, double *AVAL,
61 REAL_DH * work, int *o2n_col, Euclid_dh ctx);
62
63
64/* all non-local column indices are discarded in symbolic_row_private() */
65#undef __FUNC__
66#define __FUNC__ "iluk_mpi_bj"
67void
69{
70 START_FUNC_DH int *rp, *cval, *diag;
71 int *CVAL;
72 int i, j, len, count, col, idx = 0;
73 int *list, *marker, *fill, *tmpFill;
74 int temp, m, from = ctx->from, to = ctx->to;
75 int *n2o_row, *o2n_col;
76 int first_row, last_row;
77 double *AVAL;
78 REAL_DH *work, *aval;
79 Factor_dh F = ctx->F;
80 SubdomainGraph_dh sg = ctx->sg;
81
82 if (ctx->F == NULL)
83 {
84 SET_V_ERROR ("ctx->F is NULL");
85 }
86 if (ctx->F->rp == NULL)
87 {
88 SET_V_ERROR ("ctx->F->rp is NULL");
89 }
90
91/* printf_dh("====================== starting iluk_mpi_bj; level= %i\n\n", ctx->level);
92*/
93
94 m = F->m;
95 rp = F->rp;
96 cval = F->cval;
97 fill = F->fill;
98 diag = F->diag;
99 aval = F->aval;
100 work = ctx->work;
101
102 n2o_row = sg->n2o_row;
103 o2n_col = sg->o2n_col;
104
105 /* allocate and initialize working space */
106 list = (int *) MALLOC_DH ((m + 1) * sizeof (int));
108 marker = (int *) MALLOC_DH (m * sizeof (int));
110 tmpFill = (int *) MALLOC_DH (m * sizeof (int));
112 for (i = 0; i < m; ++i)
113 {
114 marker[i] = -1;
115 work[i] = 0.0;
116 }
117
118 /*---------- main loop ----------*/
119
120 /* global numbers of first and last locally owned rows,
121 with respect to A
122 */
123 first_row = sg->beg_row[myid_dh];
124 last_row = first_row + sg->row_count[myid_dh];
125 for (i = from; i < to; ++i)
126 {
127
128 int row = n2o_row[i]; /* local row number */
129 int globalRow = row + first_row; /* global row number */
130
131 EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
133
134 /* compute scaling value for row(i) */
135 if (ctx->isScaled)
136 {
137 compute_scaling_private (i, len, AVAL, ctx);
139 }
140
141 /* Compute symbolic factor for row(i);
142 this also performs sparsification
143 */
144 count = symbolic_row_private (i, first_row, last_row,
145 list, marker, tmpFill,
146 len, CVAL, AVAL, o2n_col, ctx);
148
149 /* Ensure adequate storage; reallocate, if necessary. */
150 if (idx + count > F->alloc)
151 {
152 Factor_dhReallocate (F, idx, count);
154 SET_INFO ("REALLOCATED from lu_mpi_bj");
155 cval = F->cval;
156 fill = F->fill;
157 aval = F->aval;
158 }
159
160 /* Copy factored symbolic row to permanent storage */
161 col = list[m];
162 while (count--)
163 {
164 cval[idx] = col;
165 fill[idx] = tmpFill[col];
166 ++idx;
167 col = list[col];
168 }
169
170 /* add row-pointer to start of next row. */
171 rp[i + 1] = idx;
172
173 /* Insert pointer to diagonal */
174 temp = rp[i];
175 while (cval[temp] != i)
176 ++temp;
177 diag[i] = temp;
178
179 /* compute numeric factor for current row */
180 numeric_row_private (i, first_row, last_row,
181 len, CVAL, AVAL, work, o2n_col, ctx);
182 CHECK_V_ERROR EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
184
185 /* Copy factored numeric row to permanent storage,
186 and re-zero work vector
187 */
188 for (j = rp[i]; j < rp[i + 1]; ++j)
189 {
190 col = cval[j];
191 aval[j] = work[col];
192 work[col] = 0.0;
193 }
194
195 /* check for zero diagonal */
196 if (!aval[diag[i]])
197 {
198 sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
200 }
201 }
202
203 FREE_DH (list);
205 FREE_DH (tmpFill);
207 FREE_DH (marker);
209
211
212
213
214/* Computes ILU(K) factor of a single row; returns fill
215 count for the row. Explicitly inserts diag if not already
216 present. On return, all column indices are local
217 (i.e, referenced to 0).
218*/
219#undef __FUNC__
220#define __FUNC__ "symbolic_row_private"
221int
222symbolic_row_private (int localRow, int beg_row, int end_row,
223 int *list, int *marker, int *tmpFill,
224 int len, int *CVAL, double *AVAL,
225 int *o2n_col, Euclid_dh ctx)
226{
227 START_FUNC_DH int level = ctx->level, m = ctx->F->m;
228 int *cval = ctx->F->cval, *diag = ctx->F->diag, *rp = ctx->F->rp;
229 int *fill = ctx->F->fill;
230 int count = 0;
231 int j, node, tmp, col, head;
232 int fill1, fill2;
233 float val;
234 double thresh = ctx->sparseTolA;
235 REAL_DH scale;
236
237 scale = ctx->scale[localRow];
238 ctx->stats[NZA_STATS] += (double) len;
239
240 /* Insert col indices in linked list, and values in work vector.
241 * List[m] points to the first (smallest) col in the linked list.
242 * Column values are adjusted from global to local numbering.
243 */
244 list[m] = m;
245 for (j = 0; j < len; ++j)
246 {
247 tmp = m;
248 col = *CVAL++;
249 val = *AVAL++;
250
251 /* throw out nonlocal columns */
252 if (col >= beg_row && col < end_row)
253 {
254 col -= beg_row; /* adjust column to local zero-based */
255 col = o2n_col[col]; /* permute column */
256 if (fabs (scale * val) > thresh || col == localRow)
257 { /* sparsification */
258 ++count;
259 while (col > list[tmp])
260 tmp = list[tmp];
261 list[col] = list[tmp];
262 list[tmp] = col;
263 tmpFill[col] = 0;
264 marker[col] = localRow;
265 }
266 }
267 }
268
269 /* insert diag if not already present */
270 if (marker[localRow] != localRow)
271 {
272/* ctx->symbolicZeroDiags += 1; */
273 tmp = m;
274 while (localRow > list[tmp])
275 tmp = list[tmp];
276 list[localRow] = list[tmp];
277 list[tmp] = localRow;
278 tmpFill[localRow] = 0;
279 marker[localRow] = localRow;
280 ++count;
281 }
282 ctx->stats[NZA_USED_STATS] += (double) count;
283
284 /* update row from previously factored rows */
285 head = m;
286 if (level > 0)
287 {
288 while (list[head] < localRow)
289 {
290 node = list[head];
291 fill1 = tmpFill[node];
292
293 if (fill1 < level)
294 {
295 for (j = diag[node] + 1; j < rp[node + 1]; ++j)
296 {
297 col = cval[j];
298 fill2 = fill1 + fill[j] + 1;
299
300 if (fill2 <= level)
301 {
302 /* if newly discovered fill entry, mark it as discovered;
303 * if entry has level <= K, add it to the linked-list.
304 */
305 if (marker[col] < localRow)
306 {
307 tmp = head;
308 marker[col] = localRow;
309 tmpFill[col] = fill2;
310 while (col > list[tmp])
311 tmp = list[tmp];
312 list[col] = list[tmp];
313 list[tmp] = col;
314 ++count; /* increment fill count */
315 }
316
317 /* if previously-discovered fill, update the entry's level. */
318 else
319 {
320 tmpFill[col] =
321 (fill2 < tmpFill[col]) ? fill2 : tmpFill[col];
322 }
323 }
324 }
325 }
326 head = list[head]; /* advance to next item in linked list */
327 }
328 }
329END_FUNC_VAL (count)}
330
331
332#undef __FUNC__
333#define __FUNC__ "numeric_row_private"
334int
335numeric_row_private (int localRow, int beg_row, int end_row,
336 int len, int *CVAL, double *AVAL,
337 REAL_DH * work, int *o2n_col, Euclid_dh ctx)
338{
339 START_FUNC_DH double pc, pv, multiplier;
340 int j, k, col, row;
341 int *rp = ctx->F->rp, *cval = ctx->F->cval;
342 int *diag = ctx->F->diag;
343 double val;
344 REAL_DH *aval = ctx->F->aval, scale;
345
346 scale = ctx->scale[localRow];
347
348 /* zero work vector */
349 /* note: indices in col[] are already permuted, and are
350 local (zero-based)
351 */
352 for (j = rp[localRow]; j < rp[localRow + 1]; ++j)
353 {
354 col = cval[j];
355 work[col] = 0.0;
356 }
357
358 /* init work vector with values from A */
359 /* (note: some values may be na due to sparsification; this is O.K.) */
360 for (j = 0; j < len; ++j)
361 {
362 col = *CVAL++;
363 val = *AVAL++;
364
365 if (col >= beg_row && col < end_row)
366 {
367 col -= beg_row; /* adjust column to local zero-based */
368 col = o2n_col[col]; /* we permute the indices from A */
369 work[col] = val * scale;
370 }
371 }
372
373 for (j = rp[localRow]; j < diag[localRow]; ++j)
374 {
375 row = cval[j];
376 pc = work[row];
377
378 if (pc != 0.0)
379 {
380 pv = aval[diag[row]];
381 multiplier = pc / pv;
382 work[row] = multiplier;
383
384 for (k = diag[row] + 1; k < rp[row + 1]; ++k)
385 {
386 col = cval[k];
387 work[col] -= (multiplier * aval[k]);
388 }
389 }
390 }
391
392 /* check for zero or too small of a pivot */
393#if 0
394 if (fabs (work[i]) <= pivotTol)
395 {
396 /* yuck! assume row scaling, and just stick in a value */
397 aval[diag[i]] = pivotFix;
398 }
399#endif
400
401END_FUNC_VAL (0)}
@ NZA_STATS
Definition: Euclid_dh.h:125
@ NZA_USED_STATS
Definition: Euclid_dh.h:127
void Factor_dhReallocate(Factor_dh F, int used, int additional)
Definition: Factor_dh.c:1185
int myid_dh
Definition: globalObjects.c:63
char msgBuf_dh[MSG_BUF_SIZE_DH]
Definition: globalObjects.c:61
#define REAL_DH
Definition: euclid_common.h:53
#define MALLOC_DH(s)
#define FREE_DH(p)
void EuclidGetRow(void *A, int row, int *len, int **ind, double **val)
Definition: getRow_dh.c:56
void EuclidRestoreRow(void *A, int row, int *len, int **ind, double **val)
Definition: getRow_dh.c:80
void compute_scaling_private(int row, int len, double *AVAL, Euclid_dh ctx)
Definition: ilu_seq.c:69
int symbolic_row_private(int localRow, int beg_row, int end_row, int *list, int *marker, int *tmpFill, int len, int *CVAL, double *AVAL, int *o2n_col, Euclid_dh ctx)
Definition: ilu_mpi_bj.c:222
void iluk_mpi_bj(Euclid_dh ctx)
Definition: ilu_mpi_bj.c:68
static int numeric_row_private(int localRow, int beg_row, int end_row, int len, int *CVAL, double *AVAL, REAL_DH *work, int *o2n_col, Euclid_dh ctx)
Definition: ilu_mpi_bj.c:335
#define SET_V_ERROR(msg)
Definition: macros_dh.h:126
#define START_FUNC_DH
Definition: macros_dh.h:181
#define CHECK_V_ERROR
Definition: macros_dh.h:138
#define SET_INFO(msg)
Definition: macros_dh.h:156
#define END_FUNC_DH
Definition: macros_dh.h:187
#define END_FUNC_VAL(retval)
Definition: macros_dh.h:208
int alloc
Definition: Factor_dh.h:74
int * fill
Definition: Factor_dh.h:72
int * cval
Definition: Factor_dh.h:70
int * diag
Definition: Factor_dh.h:73
int * rp
Definition: Factor_dh.h:69
REAL_DH * aval
Definition: Factor_dh.h:71
SubdomainGraph_dh sg
Definition: Euclid_dh.h:153
REAL_DH * scale
Definition: Euclid_dh.h:155
double stats[STATS_BINS]
Definition: Euclid_dh.h:190