| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| #include "slu_zdefs.h" |
| |
| #ifdef DEBUG |
| int num_drop_U; |
| #endif |
| |
| extern void zcopy_(int *, doublecomplex [], int *, doublecomplex [], int *); |
| |
| #if 0 |
| static doublecomplex *A; |
| static int _compare_(const void *a, const void *b) |
| { |
| register int *x = (int *)a, *y = (int *)b; |
| register double xx = z_abs1(&A[*x]), yy = z_abs1(&A[*y]); |
| if (xx > yy) return -1; |
| else if (xx < yy) return 1; |
| else return 0; |
| } |
| #endif |
| |
| int |
| ilu_zcopy_to_ucol( |
| int jcol, |
| int nseg, |
| int *segrep, |
| int *repfnz, |
| int *perm_r, |
| doublecomplex *dense, |
| int drop_rule, |
| milu_t milu, |
| double drop_tol, |
| int quota, |
| doublecomplex *sum, |
| int *nnzUj, |
| GlobalLU_t *Glu, |
| double *work |
| |
| ) |
| { |
| |
| |
| |
| int ksub, krep, ksupno; |
| int i, k, kfnz, segsze; |
| int fsupc, isub, irow; |
| int jsupno, nextu; |
| int new_next, mem_error; |
| int *xsup, *supno; |
| int *lsub, *xlsub; |
| doublecomplex *ucol; |
| int *usub, *xusub; |
| int nzumax; |
| int m; |
| register double d_max = 0.0, d_min = 1.0 / dlamch_("Safe minimum"); |
| register double tmp; |
| doublecomplex zero = {0.0, 0.0}; |
| int i_1 = 1; |
| |
| xsup = Glu->xsup; |
| supno = Glu->supno; |
| lsub = Glu->lsub; |
| xlsub = Glu->xlsub; |
| ucol = Glu->ucol; |
| usub = Glu->usub; |
| xusub = Glu->xusub; |
| nzumax = Glu->nzumax; |
| |
| *sum = zero; |
| if (drop_rule == NODROP) { |
| drop_tol = -1.0, quota = Glu->n; |
| } |
| |
| jsupno = supno[jcol]; |
| nextu = xusub[jcol]; |
| k = nseg - 1; |
| for (ksub = 0; ksub < nseg; ksub++) { |
| krep = segrep[k--]; |
| ksupno = supno[krep]; |
| |
| if ( ksupno != jsupno ) { |
| kfnz = repfnz[krep]; |
| if ( kfnz != EMPTY ) { |
| |
| fsupc = xsup[ksupno]; |
| isub = xlsub[fsupc] + kfnz - fsupc; |
| segsze = krep - kfnz + 1; |
| |
| new_next = nextu + segsze; |
| while ( new_next > nzumax ) { |
| if ((mem_error = zLUMemXpand(jcol, nextu, UCOL, &nzumax, |
| Glu)) != 0) |
| return (mem_error); |
| ucol = Glu->ucol; |
| if ((mem_error = zLUMemXpand(jcol, nextu, USUB, &nzumax, |
| Glu)) != 0) |
| return (mem_error); |
| usub = Glu->usub; |
| lsub = Glu->lsub; |
| } |
| |
| for (i = 0; i < segsze; i++) { |
| irow = lsub[isub++]; |
| tmp = z_abs1(&dense[irow]); |
| |
| |
| if (quota > 0 && tmp >= drop_tol) { |
| if (tmp > d_max) d_max = tmp; |
| if (tmp < d_min) d_min = tmp; |
| usub[nextu] = perm_r[irow]; |
| ucol[nextu] = dense[irow]; |
| nextu++; |
| } else { |
| switch (milu) { |
| case SMILU_1: |
| case SMILU_2: |
| z_add(sum, sum, &dense[irow]); |
| break; |
| case SMILU_3: |
| |
| sum->r += tmp; |
| break; |
| case SILU: |
| default: |
| break; |
| } |
| #ifdef DEBUG |
| num_drop_U++; |
| #endif |
| } |
| dense[irow] = zero; |
| } |
| |
| } |
| |
| } |
| |
| } |
| |
| xusub[jcol + 1] = nextu; |
| m = xusub[jcol + 1] - xusub[jcol]; |
| |
| |
| if (drop_rule & DROP_SECONDARY && m > quota) { |
| register double tol = d_max; |
| register int m0 = xusub[jcol] + m - 1; |
| |
| if (quota > 0) { |
| if (drop_rule & DROP_INTERP) { |
| d_max = 1.0 / d_max; d_min = 1.0 / d_min; |
| tol = 1.0 / (d_max + (d_min - d_max) * quota / m); |
| } else { |
| i_1 = xusub[jcol]; |
| for (i = 0; i < m; ++i, ++i_1) work[i] = z_abs1(&ucol[i_1]); |
| tol = dqselect(m, work, quota); |
| #if 0 |
| A = &ucol[xusub[jcol]]; |
| for (i = 0; i < m; i++) work[i] = i; |
| qsort(work, m, sizeof(int), _compare_); |
| tol = fabs(usub[xusub[jcol] + work[quota]]); |
| #endif |
| } |
| } |
| for (i = xusub[jcol]; i <= m0; ) { |
| if (z_abs1(&ucol[i]) <= tol) { |
| switch (milu) { |
| case SMILU_1: |
| case SMILU_2: |
| z_add(sum, sum, &ucol[i]); |
| break; |
| case SMILU_3: |
| sum->r += tmp; |
| break; |
| case SILU: |
| default: |
| break; |
| } |
| ucol[i] = ucol[m0]; |
| usub[i] = usub[m0]; |
| m0--; |
| m--; |
| #ifdef DEBUG |
| num_drop_U++; |
| #endif |
| xusub[jcol + 1]--; |
| continue; |
| } |
| i++; |
| } |
| } |
| |
| if (milu == SMILU_2) { |
| sum->r = z_abs1(sum); sum->i = 0.0; |
| } |
| if (milu == SMILU_3) sum->i = 0.0; |
| |
| *nnzUj += m; |
| |
| return 0; |
| } |