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"
53 static bool check_constraint_private (
Euclid_dh ctx,
int b,
int j);
55 static int symbolic_row_private (
int localRow,
56 int *list,
int *marker,
int *tmpFill,
57 int len,
int *CVAL,
double *AVAL,
60 static int numeric_row_private (
int localRow,
61 int len,
int *CVAL,
double *AVAL,
62 REAL_DH * work,
int *o2n_col,
Euclid_dh ctx,
67 #define __FUNC__ "compute_scaling_private"
69 compute_scaling_private (
int row,
int len,
double *AVAL,
Euclid_dh ctx)
71 START_FUNC_DH
double tmp = 0.0;
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]);
99 END_FUNC_VAL (max * ctxPrivate->pivotFix)}
107 #define __FUNC__ "iluk_seq"
111 START_FUNC_DH
int *rp, *cval, *diag;
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;
118 REAL_DH *work, *aval;
123 if (logFile != NULL && Parser_dhHasSwitch (parser_dh,
"-debug_ilu"))
137 SET_V_ERROR (
"subdomain graph is NULL");
140 n2o_row = ctx->sg->n2o_row;
141 o2n_col = ctx->sg->o2n_col;
142 beg_row = ctx->sg->beg_row[myid_dh];
143 beg_rowP = ctx->sg->beg_rowP[myid_dh];
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",
177 i + 1, i + 1 + sg->beg_rowP[myid_dh], ctx->level);
180 EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
186 compute_scaling_private (i, len, AVAL, ctx);
193 count = symbolic_row_private (i, list, marker, tmpFill,
194 len, CVAL, AVAL, o2n_col, ctx, debug);
198 if (idx + count > F->alloc)
200 Factor_dhReallocate (F, idx, count);
202 SET_INFO (
"REALLOCATED from ilu_seq");
213 fill[idx] = tmpFill[col];
225 while (cval[temp] != i)
233 numeric_row_private (i, len, CVAL, AVAL, work, o2n_col, ctx, debug);
234 CHECK_V_ERROR EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
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]);
251 fprintf (logFile,
"\n");
255 for (j = rp[i]; j < rp[i + 1]; ++j)
266 sprintf (msgBuf_dh,
"zero diagonal in local row %i", i + 1);
267 SET_V_ERROR (msgBuf_dh);
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"
301 START_FUNC_DH
int *rp, *cval, *diag;
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];
309 REAL_DH *work, *aval;
312 bool bj =
false, constrained =
false;
317 if (logFile != NULL && Parser_dhHasSwitch (parser_dh,
"-debug_ilu"))
323 if (!strcmp (ctx->algo_par,
"bj"))
325 constrained = !Parser_dhHasSwitch (parser_dh,
"-unconstrained");
337 n2o_row = sg->n2o_row;
338 o2n_col = sg->o2n_col;
339 row_count = sg->row_count;
341 beg_rowP = sg->beg_rowP;
342 n2o_sub = sg->n2o_sub;
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);
406 EuclidGetRow (ctx->A, row, &len, &CVAL, &AVAL);
412 compute_scaling_private (i, len, AVAL, ctx);
419 count = symbolic_row_private (i, list, marker, tmpFill,
420 len, CVAL, AVAL, o2n_col, ctx, debug);
424 if (idx + count > F->alloc)
426 Factor_dhReallocate (F, idx, count);
428 SET_INFO (
"REALLOCATED from ilu_seq");
440 if (constrained && !bj)
442 if (col >= first_row && col < end_row)
445 fill[idx] = tmpFill[col];
450 if (check_constraint_private (ctx, curBlock, 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)
500 numeric_row_private (i, len, CVAL, AVAL, work, o2n_col, ctx, debug);
501 CHECK_V_ERROR EuclidRestoreRow (ctx->A, row, &len, &CVAL, &AVAL);
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],
518 fprintf (logFile,
"\n");
524 for (j = rp[i]; j < rp[i + 1]; ++j)
535 sprintf (msgBuf_dh,
"zero diagonal in local row %i", i + 1);
536 SET_V_ERROR (msgBuf_dh);
565 #define __FUNC__ "symbolic_row_private"
567 symbolic_row_private (
int localRow,
568 int *list,
int *marker,
int *tmpFill,
569 int len,
int *CVAL,
double *AVAL,
572 START_FUNC_DH
int level = ctx->level, m = ctx->F->m;
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;
579 double thresh = ctx->sparseTolA;
582 scale = ctx->scale[localRow];
583 ctx->stats[NZA_STATS] += (double) len;
584 beg_row = ctx->sg->beg_row[myid_dh];
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;
624 ctx->stats[NZA_USED_STATS] += (double) count;
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];
677 END_FUNC_VAL (count)}
681 #define __FUNC__ "numeric_row_private"
683 numeric_row_private (
int localRow,
684 int len,
int *CVAL,
double *AVAL,
685 REAL_DH * work,
int *o2n_col,
Euclid_dh ctx,
bool debug)
687 START_FUNC_DH
double pc, pv, multiplier;
689 int *rp = ctx->F->rp, *cval = ctx->F->cval;
690 int *diag = ctx->F->diag;
693 REAL_DH *aval = ctx->F->aval, scale;
695 scale = ctx->scale[localRow];
696 beg_row = ctx->sg->beg_row[myid_dh];
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;
781 int ilut_row_private (
int localRow,
int *list,
int *o2n_col,
int *marker,
782 int len,
int *CVAL,
double *AVAL,
783 REAL_DH * work,
Euclid_dh ctx,
bool debug);
786 #define __FUNC__ "ilut_seq"
790 START_FUNC_DH
int *rp, *cval, *diag, *CVAL;
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;
796 REAL_DH *work, *aval, val;
801 if (logFile != NULL && Parser_dhHasSwitch (parser_dh,
"-debug_ilu"))
813 droptol = ctx->droptol;
817 SET_V_ERROR (
"subdomain graph is NULL");
820 n2o_row = ctx->sg->n2o_row;
821 o2n_col = ctx->sg->o2n_col;
822 beg_row = ctx->sg->beg_row[myid_dh];
823 beg_rowP = ctx->sg->beg_rowP[myid_dh];
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;
844 EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
848 compute_scaling_private (i, len, AVAL, ctx);
852 count = ilut_row_private (i, list, o2n_col, marker,
853 len, CVAL, AVAL, work, ctx, debug);
856 EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
860 if (idx + count > F->alloc)
862 Factor_dhReallocate (F, idx, count);
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);
899 SET_V_ERROR (msgBuf_dh);
906 int start = rp[from];
908 for (i = start; i < stop; ++i)
918 #define __FUNC__ "ilut_row_private"
920 ilut_row_private (
int localRow,
int *list,
int *o2n_col,
int *marker,
921 int len,
int *CVAL,
double *AVAL,
922 REAL_DH * work,
Euclid_dh ctx,
bool debug)
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;
932 double droptol = ctx->droptol;
933 double thresh = ctx->sparseTolA;
935 scale = ctx->scale[localRow];
936 ctx->stats[NZA_STATS] += (double) len;
937 beg_row = ctx->sg->beg_row[myid_dh];
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];
1018 END_FUNC_VAL (count)}
1022 #define __FUNC__ "check_constraint_private"
1024 check_constraint_private (
Euclid_dh ctx,
int p1,
int j)
1026 START_FUNC_DH
bool retval =
false;
1033 SET_ERROR (-1,
"ctx->sg == NULL");
1036 p2 = SubdomainGraph_dhFindOwner (ctx->sg, j,
true);
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)
1059 END_FUNC_VAL (retval)}