IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
ilu_mpi_pilu.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 #include "Euclid_dh.h"
44 #include "Factor_dh.h"
45 #include "Mat_dh.h"
46 #include "ilu_dh.h"
47 #include "Mem_dh.h"
48 #include "Parser_dh.h"
49 #include "Hash_dh.h"
50 #include "getRow_dh.h"
51 #include "SortedList_dh.h"
52 #include "ExternalRows_dh.h"
53 #include "SubdomainGraph_dh.h"
54 
55 
56 static void iluk_symbolic_row_private (int localRow, int len, int *CVAL,
57  double *AVAL, ExternalRows_dh extRows,
58  SortedList_dh sList, Euclid_dh ctx,
59  bool debug);
60 
61 static void iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
62  SortedList_dh slist, Euclid_dh ctx,
63  bool debug);
64 
65 #undef __FUNC__
66 #define __FUNC__ "iluk_mpi_pilu"
67 void
68 iluk_mpi_pilu (Euclid_dh ctx)
69 {
70  START_FUNC_DH int from = ctx->from, to = ctx->to;
71  int i, m;
72  int *n2o_row; /* *o2n_col; */
73  int *rp, *cval, *diag, *fill;
74  int beg_row, beg_rowP, end_rowP;
75  SubdomainGraph_dh sg = ctx->sg;
76  int *CVAL, len, idx = 0, count;
77  double *AVAL;
78  REAL_DH *aval;
79  Factor_dh F = ctx->F;
80  SortedList_dh slist = ctx->slist;
81  ExternalRows_dh extRows = ctx->extRows;
82  bool bj, noValues, debug = false;
83 
84  /* for debugging */
85  if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
86  debug = true;
87  noValues = Parser_dhHasSwitch (parser_dh, "-noValues");
88  bj = ctx->F->blockJacobi;
89 
90  m = F->m;
91  rp = F->rp;
92  cval = F->cval;
93  fill = F->fill;
94  diag = F->diag;
95  aval = F->aval;
96  /* work = ctx->work; */
97 
98  n2o_row = sg->n2o_row;
99  /* o2n_col = sg->o2n_col; */
100 
101  if (from != 0)
102  idx = rp[from];
103 
104  /* global numbers of first and last locally owned rows,
105  with respect to A
106  */
107  beg_row = sg->beg_row[myid_dh];
108  /* end_row = beg_row + sg->row_count[myid_dh]; */
109 
110  /* global number or first locally owned row, after reordering */
111  beg_rowP = sg->beg_rowP[myid_dh];
112  end_rowP = beg_rowP + sg->row_count[myid_dh];
113 
114 
115  /* loop over rows to be factored (i references local rows) */
116  for (i = from; i < to; ++i)
117  {
118 
119  int row = n2o_row[i]; /* local row number */
120  int globalRow = row + beg_row; /* global row number */
121 
122  if (debug)
123  {
124  fprintf (logFile,
125  "\nILU_pilu global: %i old_Local: %i =========================================================\n",
126  i + 1 + beg_rowP, row + 1);
127  }
128 
129  EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
130  CHECK_V_ERROR;
131 
132  if (debug)
133  {
134  int h;
135  fprintf (logFile, "ILU_pilu EuclidGetRow:\n");
136  for (h = 0; h < len; ++h)
137  fprintf (logFile, " %i %g\n", 1 + CVAL[h], AVAL[h]);
138  }
139 
140 
141  /* compute scaling value for row(i) */
142  if (ctx->isScaled)
143  {
144  compute_scaling_private (i, len, AVAL, ctx);
145  CHECK_V_ERROR;
146  }
147 
148  SortedList_dhReset (slist, i);
149  CHECK_V_ERROR;
150 
151  /* Compute symbolic factor for row(i);
152  this also performs sparsification
153  */
154  iluk_symbolic_row_private (i, len, CVAL, AVAL,
155  extRows, slist, ctx, debug);
156  CHECK_V_ERROR;
157 
158  /* enforce subdomain constraint */
159  SortedList_dhEnforceConstraint (slist, sg);
160 
161  /* compute numeric factor for row */
162  if (!noValues)
163  {
164  iluk_numeric_row_private (i, extRows, slist, ctx, debug);
165  CHECK_V_ERROR;
166  }
167 
168  EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
169  CHECK_V_ERROR;
170 
171  /* Ensure adequate storage; reallocate, if necessary. */
172  count = SortedList_dhReadCount (slist);
173  CHECK_V_ERROR;
174 
175  /* Ensure adequate storage; reallocate, if necessary. */
176  if (idx + count > F->alloc)
177  {
178  Factor_dhReallocate (F, idx, count);
179  CHECK_V_ERROR;
180  SET_INFO ("REALLOCATED from ilu_mpi_pilu");
181  cval = F->cval;
182  fill = F->fill;
183  aval = F->aval;
184  }
185 
186  /* Copy factor to permanent storage */
187  if (bj)
188  { /* for debugging: blockJacobi case */
189  int col;
190  while (count--)
191  {
192  SRecord *sr = SortedList_dhGetSmallest (slist);
193  CHECK_V_ERROR;
194  col = sr->col;
195  if (col >= beg_rowP && col < end_rowP)
196  {
197  cval[idx] = col;
198  if (noValues)
199  {
200  aval[idx] = 0.0;
201  }
202  else
203  {
204  aval[idx] = sr->val;
205  }
206  fill[idx] = sr->level;
207  ++idx;
208  }
209  }
210  }
211 
212  if (debug)
213  {
214  fprintf (logFile, "ILU_pilu ");
215  while (count--)
216  {
217  SRecord *sr = SortedList_dhGetSmallest (slist);
218  CHECK_V_ERROR;
219  cval[idx] = sr->col;
220  aval[idx] = sr->val;
221  fill[idx] = sr->level;
222  fprintf (logFile, "%i,%i,%g ; ", 1 + cval[idx], fill[idx],
223  aval[idx]);
224  ++idx;
225  }
226  fprintf (logFile, "\n");
227  }
228 
229  else
230  {
231  while (count--)
232  {
233  SRecord *sr = SortedList_dhGetSmallest (slist);
234  CHECK_V_ERROR;
235  cval[idx] = sr->col;
236  aval[idx] = sr->val;
237  fill[idx] = sr->level;
238  ++idx;
239  }
240  }
241 
242  /* add row-pointer to start of next row. */
243  rp[i + 1] = idx;
244 
245  /* Insert pointer to diagonal */
246  {
247  int temp = rp[i];
248  bool flag = true;
249  while (temp < idx)
250  {
251  if (cval[temp] == i + beg_rowP)
252  {
253  diag[i] = temp;
254  flag = false;
255  break;
256  }
257  ++temp;
258  }
259  if (flag)
260  {
261  if (logFile != NULL)
262  {
263  int k;
264  fprintf (logFile,
265  "Failed to find diag in localRow %i (globalRow %i; ct= %i)\n ",
266  1 + i, i + 1 + beg_rowP, rp[i + 1] - rp[i]);
267  for (k = rp[i]; k < rp[i + 1]; ++k)
268  {
269  fprintf (logFile, "%i ", cval[i] + 1);
270  }
271  fprintf (logFile, "\n\n");
272  }
273  sprintf (msgBuf_dh, "failed to find diagonal for localRow: %i",
274  1 + i);
275  SET_V_ERROR (msgBuf_dh);
276  }
277  }
278 /*
279  { int temp = rp[i];
280  while (cval[temp] != i+beg_row) ++temp;
281  diag[i] = temp;
282  }
283 */
284 
285  /* check for zero diagonal */
286  if (!aval[diag[i]])
287  {
288  sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
289  SET_V_ERROR (msgBuf_dh);
290  }
291 
292  }
293 
294  /* adjust to local (zero) based, if block jacobi factorization */
295  if (bj)
296  {
297  int nz = rp[m];
298  for (i = 0; i < nz; ++i)
299  cval[i] -= beg_rowP;
300  }
301 
302 END_FUNC_DH}
303 
304 
305 #undef __FUNC__
306 #define __FUNC__ "iluk_symbolic_row_private"
307 void
308 iluk_symbolic_row_private (int localRow, int len, int *CVAL,
309  double *AVAL, ExternalRows_dh extRows,
310  SortedList_dh slist, Euclid_dh ctx, bool debug)
311 {
312  START_FUNC_DH int level = ctx->level, m = ctx->m;
313  int beg_row = ctx->sg->beg_row[myid_dh];
314  int beg_rowP = ctx->sg->beg_rowP[myid_dh];
315  int *cval = ctx->F->cval, *diag = ctx->F->diag;
316  int *rp = ctx->F->rp, *fill = ctx->F->fill;
317  int j, node, col;
318  int end_rowP = beg_rowP + m;
319  int level_1, level_2;
320  int *cvalPtr, *fillPtr;
321  SRecord sr, *srPtr;
322  REAL_DH scale, *avalPtr;
323  double thresh = ctx->sparseTolA;
324  bool wasInserted;
325  int count = 0;
326 
327  scale = ctx->scale[localRow];
328  ctx->stats[NZA_STATS] += (double) len;
329 
330  /* insert col indices in sorted linked list */
331  sr.level = 0;
332  for (j = 0; j < len; ++j)
333  {
334  sr.col = CVAL[j];
335  sr.val = scale * AVAL[j];
336 /* if (fabs(sr.val) > thresh) { */
337  wasInserted = SortedList_dhPermuteAndInsert (slist, &sr, thresh);
338  CHECK_V_ERROR;
339  if (wasInserted)
340  ++count;
341 /* } */
342  if (debug)
343  {
344  fprintf (logFile, "ILU_pilu inserted from A: col= %i val= %g\n",
345  1 + CVAL[j], sr.val);
346  }
347  }
348 
349  /* ensure diagonal entry is inserted */
350  sr.val = 0.0;
351  sr.col = localRow + beg_rowP;
352  srPtr = SortedList_dhFind (slist, &sr);
353  CHECK_V_ERROR;
354  if (srPtr == NULL)
355  {
356  SortedList_dhInsert (slist, &sr);
357  CHECK_V_ERROR;
358  ++count;
359  if (debug)
360  {
361  fprintf (logFile, "ILU_pilu inserted missing diagonal: %i\n",
362  1 + localRow + beg_row);
363  }
364  }
365  ctx->stats[NZA_USED_STATS] += (double) count;
366 
367  /* update row from previously factored rows */
368  sr.val = 0.0;
369  if (level > 0)
370  {
371  while (1)
372  {
373  srPtr = SortedList_dhGetSmallestLowerTri (slist);
374  CHECK_V_ERROR;
375  if (srPtr == NULL)
376  break;
377 
378  node = srPtr->col;
379 
380  if (debug)
381  {
382  fprintf (logFile, "ILU_pilu sf updating from row: %i\n",
383  1 + srPtr->col);
384  }
385 
386  level_1 = srPtr->level;
387  if (level_1 < level)
388  {
389 
390  /* case 1: locally owned row */
391  if (node >= beg_rowP && node < end_rowP)
392  {
393  node -= beg_rowP;
394  len = rp[node + 1] - diag[node] - 1;
395  cvalPtr = cval + diag[node] + 1;
396  fillPtr = fill + diag[node] + 1;
397  }
398 
399  /* case 2: external row */
400  else
401  {
402  len = 0;
403  ExternalRows_dhGetRow (extRows, node, &len, &cvalPtr,
404  &fillPtr, &avalPtr);
405  CHECK_V_ERROR;
406  if (debug && len == 0)
407  {
408  fprintf (stderr,
409  "ILU_pilu sf failed to get extern row: %i\n",
410  1 + node);
411  }
412  }
413 
414 
415  /* merge in strict upper triangular portion of row */
416  for (j = 0; j < len; ++j)
417  {
418  col = *cvalPtr++;
419  level_2 = 1 + level_1 + *fillPtr++;
420  if (level_2 <= level)
421  {
422  /* Insert new element, or update level if already inserted. */
423  sr.col = col;
424  sr.level = level_2;
425  sr.val = 0.0;
426  SortedList_dhInsertOrUpdate (slist, &sr);
427  CHECK_V_ERROR;
428  }
429  }
430  }
431  }
432  }
433 END_FUNC_DH}
434 
435 
436 #undef __FUNC__
437 #define __FUNC__ "iluk_numeric_row_private"
438 void
439 iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
440  SortedList_dh slist, Euclid_dh ctx, bool debug)
441 {
442  START_FUNC_DH int m = ctx->m;
443  int beg_rowP = ctx->sg->beg_rowP[myid_dh];
444  int end_rowP = beg_rowP + m;
445  int len, row;
446  int *rp = ctx->F->rp, *cval = ctx->F->cval, *diag = ctx->F->diag;
447  REAL_DH *avalPtr, *aval = ctx->F->aval;
448  int *cvalPtr;
449  double multiplier, pc, pv;
450  SRecord sr, *srPtr;
451 
452  /* note: non-zero entries from A were inserted in list during iluk_symbolic_row_private */
453 
454  SortedList_dhResetGetSmallest (slist);
455  CHECK_V_ERROR;
456  while (1)
457  {
458  srPtr = SortedList_dhGetSmallestLowerTri (slist);
459  CHECK_V_ERROR;
460  if (srPtr == NULL)
461  break;
462 
463  /* update new_row's values from upper triangular portion of previously
464  factored row
465  */
466  row = srPtr->col;
467 
468  if (row >= beg_rowP && row < end_rowP)
469  {
470  int local_row = row - beg_rowP;
471 
472  len = rp[local_row + 1] - diag[local_row];
473  cvalPtr = cval + diag[local_row];
474  avalPtr = aval + diag[local_row];
475  }
476  else
477  {
478  len = 0;
479  ExternalRows_dhGetRow (extRows, row, &len, &cvalPtr,
480  NULL, &avalPtr);
481  CHECK_V_ERROR;
482  if (debug && len == 0)
483  {
484  fprintf (stderr, "ILU_pilu failed to get extern row: %i\n",
485  1 + row);
486  }
487 
488  }
489 
490  if (len)
491  {
492  /* first, form and store pivot */
493  sr.col = row;
494  srPtr = SortedList_dhFind (slist, &sr);
495  CHECK_V_ERROR;
496  if (srPtr == NULL)
497  {
498  sprintf (msgBuf_dh,
499  "find failed for sr.col = %i while factoring local row= %i \n",
500  1 + sr.col, new_row + 1);
501  SET_V_ERROR (msgBuf_dh);
502  }
503 
504  pc = srPtr->val;
505 
506  if (pc != 0.0)
507  {
508  pv = *avalPtr++;
509  --len;
510  ++cvalPtr;
511  multiplier = pc / pv;
512  srPtr->val = multiplier;
513 
514  if (debug)
515  {
516  fprintf (logFile,
517  "ILU_pilu nf updating from row: %i; multiplier = %g\n",
518  1 + srPtr->col, multiplier);
519  }
520 
521  /* second, update from strict upper triangular portion of row */
522  while (len--)
523  {
524  sr.col = *cvalPtr++;
525  sr.val = *avalPtr++;
526  srPtr = SortedList_dhFind (slist, &sr);
527  CHECK_V_ERROR;
528  if (srPtr != NULL)
529  {
530  srPtr->val -= (multiplier * sr.val);
531  }
532  }
533  }
534 
535  else
536  {
537  if (debug)
538  {
539  fprintf (logFile,
540  "ILU_pilu NO UPDATE from row: %i; srPtr->val = 0.0\n",
541  1 + srPtr->col);
542  }
543  }
544 
545  }
546  }
547 END_FUNC_DH}