IFPACK Development
Loading...
Searching...
No Matches
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
56static 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
61static 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"
67void
68iluk_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
302END_FUNC_DH}
303
304
305#undef __FUNC__
306#define __FUNC__ "iluk_symbolic_row_private"
307void
308iluk_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 }
433END_FUNC_DH}
434
435
436#undef __FUNC__
437#define __FUNC__ "iluk_numeric_row_private"
438void
439iluk_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 }
547END_FUNC_DH}