43 #include "Euclid_dh.h"
44 #include "Factor_dh.h"
48 #include "Parser_dh.h"
50 #include "getRow_dh.h"
51 #include "SortedList_dh.h"
52 #include "ExternalRows_dh.h"
53 #include "SubdomainGraph_dh.h"
56 static void iluk_symbolic_row_private (
int localRow,
int len,
int *CVAL,
61 static void iluk_numeric_row_private (
int new_row,
ExternalRows_dh extRows,
66 #define __FUNC__ "iluk_mpi_pilu"
70 START_FUNC_DH
int from = ctx->from, to = ctx->to;
73 int *rp, *cval, *diag, *fill;
74 int beg_row, beg_rowP, end_rowP;
76 int *CVAL, len, idx = 0, count;
82 bool bj, noValues, debug =
false;
85 if (logFile != NULL && Parser_dhHasSwitch (parser_dh,
"-debug_ilu"))
87 noValues = Parser_dhHasSwitch (parser_dh,
"-noValues");
88 bj = ctx->F->blockJacobi;
98 n2o_row = sg->n2o_row;
107 beg_row = sg->beg_row[myid_dh];
111 beg_rowP = sg->beg_rowP[myid_dh];
112 end_rowP = beg_rowP + sg->row_count[myid_dh];
116 for (i = from; i < to; ++i)
119 int row = n2o_row[i];
120 int globalRow = row + beg_row;
125 "\nILU_pilu global: %i old_Local: %i =========================================================\n",
126 i + 1 + beg_rowP, row + 1);
129 EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
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]);
144 compute_scaling_private (i, len, AVAL, ctx);
148 SortedList_dhReset (slist, i);
154 iluk_symbolic_row_private (i, len, CVAL, AVAL,
155 extRows, slist, ctx, debug);
159 SortedList_dhEnforceConstraint (slist, sg);
164 iluk_numeric_row_private (i, extRows, slist, ctx, debug);
168 EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
172 count = SortedList_dhReadCount (slist);
176 if (idx + count > F->alloc)
178 Factor_dhReallocate (F, idx, count);
180 SET_INFO (
"REALLOCATED from ilu_mpi_pilu");
192 SRecord *sr = SortedList_dhGetSmallest (slist);
195 if (col >= beg_rowP && col < end_rowP)
206 fill[idx] = sr->level;
214 fprintf (logFile,
"ILU_pilu ");
217 SRecord *sr = SortedList_dhGetSmallest (slist);
221 fill[idx] = sr->level;
222 fprintf (logFile,
"%i,%i,%g ; ", 1 + cval[idx], fill[idx],
226 fprintf (logFile,
"\n");
233 SRecord *sr = SortedList_dhGetSmallest (slist);
237 fill[idx] = sr->level;
251 if (cval[temp] == i + beg_rowP)
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)
269 fprintf (logFile,
"%i ", cval[i] + 1);
271 fprintf (logFile,
"\n\n");
273 sprintf (msgBuf_dh,
"failed to find diagonal for localRow: %i",
275 SET_V_ERROR (msgBuf_dh);
288 sprintf (msgBuf_dh,
"zero diagonal in local row %i", i + 1);
289 SET_V_ERROR (msgBuf_dh);
298 for (i = 0; i < nz; ++i)
306 #define __FUNC__ "iluk_symbolic_row_private"
308 iluk_symbolic_row_private (
int localRow,
int len,
int *CVAL,
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;
318 int end_rowP = beg_rowP + m;
319 int level_1, level_2;
320 int *cvalPtr, *fillPtr;
322 REAL_DH scale, *avalPtr;
323 double thresh = ctx->sparseTolA;
327 scale = ctx->scale[localRow];
328 ctx->stats[NZA_STATS] += (double) len;
332 for (j = 0; j < len; ++j)
335 sr.val = scale * AVAL[j];
337 wasInserted = SortedList_dhPermuteAndInsert (slist, &sr, thresh);
344 fprintf (logFile,
"ILU_pilu inserted from A: col= %i val= %g\n",
345 1 + CVAL[j], sr.val);
351 sr.col = localRow + beg_rowP;
352 srPtr = SortedList_dhFind (slist, &sr);
356 SortedList_dhInsert (slist, &sr);
361 fprintf (logFile,
"ILU_pilu inserted missing diagonal: %i\n",
362 1 + localRow + beg_row);
365 ctx->stats[NZA_USED_STATS] += (double) count;
373 srPtr = SortedList_dhGetSmallestLowerTri (slist);
382 fprintf (logFile,
"ILU_pilu sf updating from row: %i\n",
386 level_1 = srPtr->level;
391 if (node >= beg_rowP && node < end_rowP)
394 len = rp[node + 1] - diag[node] - 1;
395 cvalPtr = cval + diag[node] + 1;
396 fillPtr = fill + diag[node] + 1;
403 ExternalRows_dhGetRow (extRows, node, &len, &cvalPtr,
406 if (debug && len == 0)
409 "ILU_pilu sf failed to get extern row: %i\n",
416 for (j = 0; j < len; ++j)
419 level_2 = 1 + level_1 + *fillPtr++;
420 if (level_2 <= level)
426 SortedList_dhInsertOrUpdate (slist, &sr);
437 #define __FUNC__ "iluk_numeric_row_private"
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;
446 int *rp = ctx->F->rp, *cval = ctx->F->cval, *diag = ctx->F->diag;
447 REAL_DH *avalPtr, *aval = ctx->F->aval;
449 double multiplier, pc, pv;
454 SortedList_dhResetGetSmallest (slist);
458 srPtr = SortedList_dhGetSmallestLowerTri (slist);
468 if (row >= beg_rowP && row < end_rowP)
470 int local_row = row - beg_rowP;
472 len = rp[local_row + 1] - diag[local_row];
473 cvalPtr = cval + diag[local_row];
474 avalPtr = aval + diag[local_row];
479 ExternalRows_dhGetRow (extRows, row, &len, &cvalPtr,
482 if (debug && len == 0)
484 fprintf (stderr,
"ILU_pilu failed to get extern row: %i\n",
494 srPtr = SortedList_dhFind (slist, &sr);
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);
511 multiplier = pc / pv;
512 srPtr->val = multiplier;
517 "ILU_pilu nf updating from row: %i; multiplier = %g\n",
518 1 + srPtr->col, multiplier);
526 srPtr = SortedList_dhFind (slist, &sr);
530 srPtr->val -= (multiplier * sr.val);
540 "ILU_pilu NO UPDATE from row: %i; srPtr->val = 0.0\n",