56 int *list,
int *marker,
int *tmpFill,
57 int len,
int *CVAL,
double *AVAL,
61 int len,
int *CVAL,
double *AVAL,
67 #define __FUNC__ "compute_scaling_private"
74 for (j = 0; j < len; ++j)
75 tmp =
MAX (tmp, fabs (AVAL[j]));
78 ctx->
scale[row] = 1.0 / tmp;
86 #define __FUNC__ "fixPivot_private"
88 fixPivot_private (
int row,
int len,
float *vals)
94 for (i = 0; i < len; ++i)
96 float tmp = fabs (vals[i]);
107 #define __FUNC__ "iluk_seq"
113 int i, j, len, count, col, idx = 0;
114 int *list, *marker, *fill, *tmpFill;
115 int temp, m, from = ctx->
from, to = ctx->
to;
116 int *n2o_row, *o2n_col, beg_row, beg_rowP;
146 list = (
int *)
MALLOC_DH ((m + 1) *
sizeof (int));
148 marker = (
int *)
MALLOC_DH (m *
sizeof (
int));
150 tmpFill = (
int *)
MALLOC_DH (m *
sizeof (
int));
152 for (i = 0; i < m; ++i)
156 for (i = 0; i < m; ++i)
165 for (i = from; i < to; ++i)
167 int row = n2o_row[i];
168 int globalRow = row + beg_row;
176 "ILU_seq ================================= starting local row: %i, (global= %i) level= %i\n",
194 len, CVAL, AVAL, o2n_col, ctx, debug);
198 if (idx + count > F->
alloc)
202 SET_INFO (
"REALLOCATED from ilu_seq");
213 fill[idx] = tmpFill[col];
225 while (cval[temp] != i)
242 fprintf (
logFile,
"ILU_seq: ");
243 for (j = rp[i]; j < rp[i + 1]; ++j)
248 fprintf (
logFile,
"%i,%i,%g ; ", 1 + cval[j], fill[j], aval[j]);
255 for (j = rp[i]; j < rp[i + 1]; ++j)
266 sprintf (
msgBuf_dh,
"zero diagonal in local row %i", i + 1);
281 int start = rp[from];
283 for (i = start; i < stop; ++i)
290 for (i = to + 1; i < m; ++i)
297 #define __FUNC__ "iluk_seq_block"
303 int h, i, j, len, count, col, idx = 0;
304 int *list, *marker, *fill, *tmpFill;
306 int *n2o_row, *o2n_col, *beg_rowP, *n2o_sub, blocks;
307 int *row_count, *dummy = NULL, dummy2[1];
312 bool bj =
false, constrained =
false;
348 dummy = (
int *)
MALLOC_DH (m *
sizeof (
int));
350 for (i = 0; i < m; ++i)
363 list = (
int *)
MALLOC_DH ((m + 1) *
sizeof (int));
365 marker = (
int *)
MALLOC_DH (m *
sizeof (
int));
367 tmpFill = (
int *)
MALLOC_DH (m *
sizeof (
int));
369 for (i = 0; i < m; ++i)
373 for (i = 0; i < m; ++i)
378 for (h = 0; h < blocks; ++h)
381 int curBlock = n2o_sub[h];
382 int first_row = beg_rowP[curBlock];
383 int end_row = first_row + row_count[curBlock];
387 fprintf (
logFile,
"\n\nILU_seq BLOCK: %i @@@@@@@@@@@@@@@ \n",
391 for (i = first_row; i < end_row; ++i)
393 int row = n2o_row[i];
399 "ILU_seq global: %i local: %i =================================\n",
400 1 + gr, 1 + i - first_row);
420 len, CVAL, AVAL, o2n_col, ctx, debug);
424 if (idx + count > F->
alloc)
428 SET_INFO (
"REALLOCATED from ilu_seq");
440 if (constrained && !bj)
442 if (col >= first_row && col < end_row)
445 fill[idx] = tmpFill[col];
453 fill[idx] = tmpFill[col];
467 if (col >= first_row && col < end_row)
470 fill[idx] = tmpFill[col];
484 fill[idx] = tmpFill[col];
495 while (cval[temp] != i)
509 fprintf (
logFile,
"ILU_seq: ");
510 for (j = rp[i]; j < rp[i + 1]; ++j)
515 fprintf (
logFile,
"%i,%i,%g ; ", 1 + cval[j], fill[j],
524 for (j = rp[i]; j < rp[i + 1]; ++j)
535 sprintf (
msgBuf_dh,
"zero diagonal in local row %i", i + 1);
565 #define __FUNC__ "symbolic_row_private"
568 int *list,
int *marker,
int *tmpFill,
569 int len,
int *CVAL,
double *AVAL,
573 int *cval = ctx->
F->
cval, *diag = ctx->
F->
diag, *rp = ctx->
F->
rp;
574 int *fill = ctx->
F->
fill;
576 int j, node, tmp, col, head;
577 int fill1, fill2, beg_row;
582 scale = ctx->
scale[localRow];
591 for (j = 0; j < len; ++j)
600 if (fabs (val) > thresh || col == localRow)
603 while (col > list[tmp])
605 list[col] = list[tmp];
608 marker[col] = localRow;
613 if (marker[localRow] != localRow)
616 while (localRow > list[tmp])
618 list[localRow] = list[tmp];
619 list[tmp] = localRow;
620 tmpFill[localRow] = 0;
621 marker[localRow] = localRow;
630 while (list[head] < localRow)
633 fill1 = tmpFill[node];
637 fprintf (
logFile,
"ILU_seq sf updating from row: %i\n",
643 for (j = diag[node] + 1; j < rp[node + 1]; ++j)
646 fill2 = fill1 + fill[j] + 1;
653 if (marker[col] < localRow)
656 marker[col] = localRow;
657 tmpFill[col] = fill2;
658 while (col > list[tmp])
660 list[col] = list[tmp];
669 (fill2 < tmpFill[col]) ? fill2 : tmpFill[col];
681 #define __FUNC__ "numeric_row_private"
684 int len,
int *CVAL,
double *AVAL,
689 int *rp = ctx->
F->
rp, *cval = ctx->
F->
cval;
690 int *diag = ctx->
F->
diag;
695 scale = ctx->
scale[localRow];
700 for (j = rp[localRow]; j < rp[localRow + 1]; ++j)
708 for (j = 0; j < len; ++j)
714 work[col] = val * scale;
723 for (j = rp[localRow]; j < diag[localRow]; ++j)
729 pv = aval[diag[row]];
737 if (pc != 0.0 && pv != 0.0)
739 multiplier = pc / pv;
740 work[row] = multiplier;
745 "ILU_seq nf updating from row: %i; multiplier= %g\n",
746 1 + row, multiplier);
749 for (k = diag[row] + 1; k < rp[row + 1]; ++k)
752 work[col] -= (multiplier * aval[k]);
760 "ILU_seq nf NO UPDATE from row %i; pc = %g; pv = %g\n",
768 if (fabs (work[i]) <= pivotTol)
771 aval[diag[i]] = pivotFix;
782 int len,
int *CVAL,
double *AVAL,
786 #define __FUNC__ "ilut_seq"
791 int i, len, count, col, idx = 0;
793 int temp, m, from, to;
794 int *n2o_row, *o2n_col, beg_row, beg_rowP;
795 double *AVAL, droptol;
827 list = (
int *)
MALLOC_DH ((m + 1) *
sizeof (int));
829 marker = (
int *)
MALLOC_DH (m *
sizeof (
int));
831 for (i = 0; i < m; ++i)
836 for (i = 0; i < m; ++i)
840 for (i = from; i < to; ++i)
842 int row = n2o_row[i];
843 int globalRow = row + beg_row;
853 len, CVAL, AVAL, work, ctx, debug);
860 if (idx + count > F->
alloc)
864 SET_INFO (
"REALLOCATED from ilu_seq");
877 if (col == i || fabs (val) > droptol)
891 while (cval[temp] != i)
898 sprintf (
msgBuf_dh,
"zero diagonal in local row %i", i + 1);
906 int start = rp[from];
908 for (i = start; i < stop; ++i)
918 #define __FUNC__ "ilut_row_private"
921 int len,
int *CVAL,
double *AVAL,
925 int j, col, m = ctx->
m, *rp = F->rp, *cval = F->cval;
926 int tmp, *diag = F->diag;
928 int count = 0, beg_row;
930 double mult, *aval = F->aval;
931 double scale, pv, pc;
935 scale = ctx->
scale[localRow];
945 for (j = 0; j < len; ++j)
954 if (fabs (val) > thresh || col == localRow)
957 while (col > list[tmp])
959 list[col] = list[tmp];
962 marker[col] = localRow;
967 if (marker[localRow] != localRow)
970 while (localRow > list[tmp])
972 list[localRow] = list[tmp];
973 list[tmp] = localRow;
974 marker[localRow] = localRow;
980 while (list[head] < localRow)
982 int row = list[head];
988 pv = aval[diag[row]];
992 if (fabs (mult) > droptol)
996 for (j = diag[row] + 1; j < rp[row + 1]; ++j)
999 work[col] -= (mult * aval[j]);
1002 if (marker[col] < localRow)
1004 marker[col] = localRow;
1006 while (col > list[tmp])
1008 list[col] = list[tmp];
1022 #define __FUNC__ "check_constraint_private"
1039 nabors = sg->
adj + sg->
ptrs[p1];
1040 count = sg->
ptrs[p1 + 1] - sg->
ptrs[p1];
1048 for (i = 0; i < count; ++i)
1052 if (nabors[i] == p2)
void ilut_seq(Euclid_dh ctx)
void EuclidGetRow(void *A, int row, int *len, int **ind, double **val)
void compute_scaling_private(int row, int len, double *AVAL, Euclid_dh ctx)
#define END_FUNC_VAL(retval)
char algo_par[MAX_OPT_LEN]
#define SET_ERROR(retval, msg)
bool Parser_dhHasSwitch(Parser_dh p, char *s)
static bool check_constraint_private(Euclid_dh ctx, int b, int j)
static int numeric_row_private(int localRow, int len, int *CVAL, double *AVAL, REAL_DH *work, int *o2n_col, Euclid_dh ctx, bool debug)
int ilut_row_private(int localRow, int *list, int *o2n_col, int *marker, int len, int *CVAL, double *AVAL, REAL_DH *work, Euclid_dh ctx, bool debug)
void iluk_seq_block(Euclid_dh ctx)
int SubdomainGraph_dhFindOwner(SubdomainGraph_dh s, int idx, bool permuted)
void iluk_seq(Euclid_dh ctx)
int symbolic_row_private(int localRow, int beg_row, int end_row, int *list, int *marker, int *tmpFill, int len, int *CVAL, double *AVAL, int *o2n_col, Euclid_dh ctx)
void EuclidRestoreRow(void *A, int row, int *len, int **ind, double **val)
void Factor_dhReallocate(Factor_dh F, int used, int additional)
char msgBuf_dh[MSG_BUF_SIZE_DH]