IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
ilu_mpi_bj.c
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 
54 int 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 
59 static 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"
67 void
68 iluk_mpi_bj (Euclid_dh ctx)
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));
107  CHECK_V_ERROR;
108  marker = (int *) MALLOC_DH (m * sizeof (int));
109  CHECK_V_ERROR;
110  tmpFill = (int *) MALLOC_DH (m * sizeof (int));
111  CHECK_V_ERROR;
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);
132  CHECK_V_ERROR;
133 
134  /* compute scaling value for row(i) */
135  if (ctx->isScaled)
136  {
137  compute_scaling_private (i, len, AVAL, ctx);
138  CHECK_V_ERROR;
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);
147  CHECK_V_ERROR;
148 
149  /* Ensure adequate storage; reallocate, if necessary. */
150  if (idx + count > F->alloc)
151  {
152  Factor_dhReallocate (F, idx, count);
153  CHECK_V_ERROR;
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);
183  CHECK_V_ERROR;
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);
199  SET_V_ERROR (msgBuf_dh);
200  }
201  }
202 
203  FREE_DH (list);
204  CHECK_V_ERROR;
205  FREE_DH (tmpFill);
206  CHECK_V_ERROR;
207  FREE_DH (marker);
208  CHECK_V_ERROR;
209 
210 END_FUNC_DH}
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"
221 int
222 symbolic_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  }
329 END_FUNC_VAL (count)}
330 
331 
332 #undef __FUNC__
333 #define __FUNC__ "numeric_row_private"
334 int
335 numeric_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 
401 END_FUNC_VAL (0)}