49 #define MIN(a,b) ((a)<=(b) ? (a) : (b))
50 #define MAX(a,b) ((a)>=(b) ? (a) : (b))
51 #define ABS(a) ((a)>=0 ? (a) : -(a))
53 #define SHORTCUT(p, a, ja, ia) \
65 a->
val =
new double[nnz];
66 a->
col =
new int[nnz];
67 a->
ptr =
new int[n+1];
80 static void qsplit(
double *a,
int *ind,
int n,
int ncut)
83 int itmp, first, last, mid;
89 if (ncut < first || ncut > last)
97 for (j=first+1; j<=last; j++)
99 if (
ABS(a[j]) > abskey)
117 ind[mid] = ind[first];
137 const int *ia,
const int *ja,
const double *a,
141 const double *multipliers,
148 int j, i, isj, iej, row;
169 mult = multipliers[ifirst2[j]] * d[j];
174 for (i=isj; i<=iej; i++)
182 if ((pos = marker[row]) != -1)
185 ta[pos] -= mult*a[i];
191 ta[talen] = -mult*a[i];
192 marker[row] = talen++;
212 int j, isj, iej, iptr;
222 if (isj != 0 && isj <= iej)
230 list[iptr] = list[ja[isj]];
231 list[ja[isj]] = iptr;
248 list[k] = list[iptr];
280 int (*getcol)(
void * A,
int col,
int ** nentries,
double * val,
int * ind),
281 int (*getdiag)(
void *A,
double * diag),
296 double *work_l =
new double[n];
297 int *ind_l =
new int[n];
301 int *list_l =
new int[n];
302 int *ifirst_l =
new int[n];
305 int *marker_l = ifirst_l;
308 double *al;
int *jal, *ial;
309 double *l;
int *jl, *il;
313 double *diag =
new double[n];
354 getcol(A, k, len_l, work_l, ind_l);
355 for (j=0; j<len_l; j++)
357 norm_l +=
ABS(work_l[j]);
358 marker_l[ind_l[j]] = j;
362 for (j=ial[k]; j<ial[k+1]; j++)
365 work_l[len_l] = al[j];
366 norm_l +=
ABS(al[j]);
367 ind_l[len_l] = index;
368 marker_l[index] = len_l++;
371 norm_l = (len_l == 0) ? 0.0 : norm_l/((
double) len_l);
376 diag, marker_l, work_l, ind_l, &len_l);
378 for (j=0; j<len_l; j++)
379 work_l[j] /= diag[k];
382 for (j=0; j<len_l; j++)
384 if (
ABS(work_l[j]) < droptol * norm_l)
387 marker_l[ind_l[j]] = -1;
391 work_l[i] = work_l[j];
405 count_l =
MIN(len_l, lfil);
406 qsplit(work_l, ind_l, len_l, count_l);
409 for (j=0; j<count_l; j++)
426 for (j=0; j<len_l; j++)
427 marker_l[ind_l[j]] = -1;
433 for (j=0; j<count_l; j++)
436 diag[index] -= work_l[j] * work_l[j] * diag[k];
#define SHORTCUT(p, a, ja, ia)
static void update_lists(int k, const int *ia, const int *ja, int *ifirst, int *list)
void Ifpack_AIJMatrix_alloc(Ifpack_AIJMatrix *a, int n, int nnz)
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
static void update_lists_newcol(int k, int isk, int iptr, int *ifirst, int *list)
static void qsplit(double *a, int *ind, int n, int ncut)
void ifpack_multilist_sort(int *const pbase, double *const daux, int total_elems)
static void update_column(int k, const int *ia, const int *ja, const double *a, const int *ifirst, const int *ifirst2, const int *list2, const double *multipliers, const double *d, int *marker, double *ta, int *itcol, int *ptalen)
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)