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"
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,
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);
66 #define __FUNC__ "iluk_mpi_bj"
70 START_FUNC_DH
int *rp, *cval, *diag;
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;
84 SET_V_ERROR (
"ctx->F is NULL");
86 if (ctx->F->rp == NULL)
88 SET_V_ERROR (
"ctx->F->rp is NULL");
102 n2o_row = sg->n2o_row;
103 o2n_col = sg->o2n_col;
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)
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)
128 int row = n2o_row[i];
129 int globalRow = row + first_row;
131 EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
137 compute_scaling_private (i, len, AVAL, ctx);
144 count = symbolic_row_private (i, first_row, last_row,
145 list, marker, tmpFill,
146 len, CVAL, AVAL, o2n_col, ctx);
150 if (idx + count > F->alloc)
152 Factor_dhReallocate (F, idx, count);
154 SET_INFO (
"REALLOCATED from lu_mpi_bj");
165 fill[idx] = tmpFill[col];
175 while (cval[temp] != i)
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);
188 for (j = rp[i]; j < rp[i + 1]; ++j)
198 sprintf (msgBuf_dh,
"zero diagonal in local row %i", i + 1);
199 SET_V_ERROR (msgBuf_dh);
220 #define __FUNC__ "symbolic_row_private"
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,
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;
231 int j, node, tmp, col, head;
234 double thresh = ctx->sparseTolA;
237 scale = ctx->scale[localRow];
238 ctx->stats[NZA_STATS] += (double) len;
245 for (j = 0; j < len; ++j)
252 if (col >= beg_row && col < end_row)
256 if (fabs (scale * val) > thresh || col == localRow)
259 while (col > list[tmp])
261 list[col] = list[tmp];
264 marker[col] = localRow;
270 if (marker[localRow] != localRow)
274 while (localRow > list[tmp])
276 list[localRow] = list[tmp];
277 list[tmp] = localRow;
278 tmpFill[localRow] = 0;
279 marker[localRow] = localRow;
282 ctx->stats[NZA_USED_STATS] += (double) count;
288 while (list[head] < localRow)
291 fill1 = tmpFill[node];
295 for (j = diag[node] + 1; j < rp[node + 1]; ++j)
298 fill2 = fill1 + fill[j] + 1;
305 if (marker[col] < localRow)
308 marker[col] = localRow;
309 tmpFill[col] = fill2;
310 while (col > list[tmp])
312 list[col] = list[tmp];
321 (fill2 < tmpFill[col]) ? fill2 : tmpFill[col];
329 END_FUNC_VAL (count)}
333 #define __FUNC__ "numeric_row_private"
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)
339 START_FUNC_DH
double pc, pv, multiplier;
341 int *rp = ctx->F->rp, *cval = ctx->F->cval;
342 int *diag = ctx->F->diag;
344 REAL_DH *aval = ctx->F->aval, scale;
346 scale = ctx->scale[localRow];
352 for (j = rp[localRow]; j < rp[localRow + 1]; ++j)
360 for (j = 0; j < len; ++j)
365 if (col >= beg_row && col < end_row)
369 work[col] = val * scale;
373 for (j = rp[localRow]; j < diag[localRow]; ++j)
380 pv = aval[diag[row]];
381 multiplier = pc / pv;
382 work[row] = multiplier;
384 for (k = diag[row] + 1; k < rp[row + 1]; ++k)
387 work[col] -= (multiplier * aval[k]);
394 if (fabs (work[i]) <= pivotTol)
397 aval[diag[i]] = pivotFix;