63 for(i = 1; i < *size; ++i)
64 for(j = 0; j < i; ++j)
65 assert(heap[i] != heap[j]);
77 e = heap[s = --(*size)];
80 for(j = 0, i = 1; i < s; i = 2 * j + 1)
113 if(i < *size && e < heap[i])
148 for(i = 1; i < *size; ++i)
149 for(j = 0; j < i; ++j)
150 assert(heap[i] != heap[j]);
162 e = heap[s = --(*size)];
165 for(j = 0, i = 1; i < s; i = 2 * j + 1)
198 if(i < *size && e > heap[i])
260 for(
int i = 0; i <
thedim; ++i)
278 diag[p_row] = 1.0 / val;
302 for(ring = list->
next; ring != list; ring = ring->
next)
306 if(l_rbeg[l_row] != n)
312 assert(l_rlen[l_row] <= l_rmax[l_row]);
314 l_rmax[l_row] = l_rlen[l_row];
315 j = i + l_rlen[l_row];
317 for(; i < j; ++i, ++n)
320 l_ridx[n] = l_ridx[i];
321 l_rval[n] = l_rval[i];
328 goto terminatePackRows;
333 l_rmax[l_row] = l_rlen[l_row];
360 for(ring = list->
next; ring != list; ring = ring->
next)
371 cmax[colno] = clen[colno];
384 goto terminatePackColumns;
389 cmax[colno] = clen[colno];
392 terminatePackColumns :
403 assert(
u.
row.
max[p_row] < len);
407 int delta = len -
u.
row.
max[p_row];
412 delta = len -
u.
row.
max[p_row];
422 &&
"ERROR: could not allocate memory for row file");
445 &&
"ERROR: could not allocate memory for row file");
460 for(; i < k; ++i, ++j)
489 for(ring = list->
next; ring != list; ring = ring->
next)
493 if(l_cbeg[l_col] != n)
500 l_cmax[l_col] = l_clen[l_col];
501 j = i + l_clen[l_col];
504 l_cidx[n++] = l_cidx[i];
510 goto terminatePackColumns;
515 l_cmax[l_col] = l_clen[l_col];
518 terminatePackColumns :
529 assert(
u.
col.
max[p_col] < len);
533 int delta = len -
u.
col.
max[p_col];
538 delta = len -
u.
col.
max[p_col];
548 &&
"ERROR: could not allocate memory for column file");
572 &&
"ERROR: could not allocate memory for column file");
598 assert(
u.
col.
max[p_col] < len);
602 int delta = len -
u.
col.
max[p_col];
607 delta = len -
u.
col.
max[p_col];
615 &&
"ERROR: could not allocate memory for column file");
638 &&
"ERROR: could not allocate memory for column file");
694 int i, j, k, h, m, n;
726 for(i += j - 1; i >= j; --i)
732 while(ridx[k] != p_col)
757 if(num > cmax[p_col])
768 for(j = 0; j < num; ++j)
780 assert(k - cbeg[p_col] < cmax[p_col]);
787 if(rmax[i] <= rlen[i])
794 h = rbeg[i] + (rlen[i])++;
806 nzCnt += (clen[p_col] = k - cbeg[p_col]);
823 for(i = 0; i < dim; ++i)
836 clen[p_col] = k - cbeg[p_col];
845 assert(k - cbeg[p_col] < cmax[p_col]);
852 if(rmax[i] <= rlen[i])
859 h = rbeg[i] + (rlen[i])++;
871 nzCnt += (clen[p_col] = k - cbeg[p_col]);
873 if(cbeg[p_col] + cmax[p_col] ==
u.
col.
used)
876 cmax[p_col] = clen[p_col];
892 memmove(&rorig[c], &rorig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
896 for(i = c; i <= r; ++i)
904 memmove(&corig[c], &corig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
908 for(i = c; i <= r; ++i)
920 if(i < verySparseFactor * (dim - c))
937 for(i += j - 1; i >= j; --i)
942 m = --(clen[k]) + cbeg[k];
944 for(h = m; cidx[h] != rowno; --h)
960 assert((num == 0) || (nonz != 0));
968 for(i = 0; i < num; ++i)
969 assert(p_work[corig[nonz[i]]] != 0);
979 assert(p_work[k] != 0);
983 x = p_work[k] *
diag[n];
1032 diag[rowno] = 1 / x;
1039 if(rmax[rowno] < num)
1053 for(i = 0; i < num; ++i)
1077 if(clen[j] >= cmax[j])
1084 cval[cbeg[j] + clen[j]] = x;
1086 cidx[cbeg[j] + clen[j]++] = rowno;
1090 rlen[rowno] = n - rbeg[rowno];
1096 for(i += j - 1; i >= j; --i)
1099 p_work[k] = rval[i];
1100 m = --(clen[k]) + cbeg[k];
1102 for(h = m; cidx[h] != rowno; --h)
1118 for(i = c; i < r; ++i)
1125 x = p_work[k] *
diag[n];
1139 p_work[ridx[j]] -= x * rval[j];
1162 diag[rowno] = 1 / x;
1172 for(i = r + 1; i < dim; ++i)
1173 if(p_work[corig[i]] != 0)
1190 for(i = r + 1; i < dim; ++i)
1208 if(clen[j] >= cmax[j])
1215 cval[cbeg[j] + clen[j]] = x;
1217 cidx[cbeg[j] + clen[j]++] = rowno;
1221 rlen[rowno] = n - rbeg[rowno];
1231 i = rbeg[rowno] + --(rlen[rowno]);
1232 diag[rowno] = 1 / rval[i];
1234 for(j = i = --(clen[p_col]) + cbeg[p_col]; cidx[i] != rowno; --i)
1259 assert(p_work[p_col] != 0);
1260 rezi = 1 / p_work[p_col];
1268 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1271 lval[ll] = rezi * p_work[j];
1278 lval[ll] = 1 - rezi;
1281 for(--i; i >= 0; --i)
1285 lval[ll] = x = rezi * p_work[j];
1305 assert(p_work[p_col] != 0);
1306 rezi = 1 / p_work[p_col];
1312 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1315 lval[ll] = rezi * p_work[j];
1321 lval[ll] = 1 - rezi;
1324 for(--i; i >= 0; --i)
1328 lval[ll] = x = rezi * p_work[j];
1370 Dring* rring, *lastrring;
1371 Dring* cring, *lastcring;
1381 for(
int i = 0; i <
thedim; i++)
1386 for(
int i = 0; i < thedim; i++)
1397 for(
int j = 0; j < k; ++j)
1429 lastrring->
next = rring;
1437 lastcring->
next = cring;
1441 for(
int i = 0; i <
thedim; i++)
1447 rring->
prev = lastrring;
1448 lastrring->
next = rring;
1453 cring->
prev = lastcring;
1454 lastcring->
next = cring;
1478 for(
int i = 0; i <
thedim; i++)
1488 for(
int j = 0; j < psv->
size() && nnonzeros <= 1; j++)
1490 if(psv->
value(j) != 0)
1502 else if(nnonzeros == 1)
1508 for(j = 0; psv->
value(j) == 0; j++)
1511 assert(j < psv->size());
1544 assert(nnonzeros >= 2);
1547 for(
int j = 0; j < psv->
size(); j++)
1573 assert(nnonzeros >= 2);
1594 int p_col, p_row, newrow;
1614 for(j = 0; j < len; ++j)
1623 for(k = n;
u.
col.
idx[k] != p_row; ++k)
1641 if(rperm[newrow] >= 0)
1651 for(k = n;
u.
row.
idx[k] != p_col; --k)
1686 int p_row, p_col, len, rs, lk;
1695 for(i = 0; i <
thedim; ++i)
1697 if(rperm[i] < 0 &&
u.
row.
len[i] == 1)
1736 for(j = k;
u.
row.
idx[j] != p_col; --j)
1795 for(i = 0; i <
thedim; ++i)
1880 for(; (r = idx[i]) != prow; ++i)
1887 for(j = k;
u.
row.
idx[j] != pcol; --j)
1918 assert(i < len &&
"ERROR: pivot column does not contain pivot row");
1920 for(++i; i < len; ++i)
1928 for(j = k;
u.
row.
idx[j] != pcol; --j)
1990 for(i = j; (c =
u.
row.
idx[i]) != pcol; --i)
1994 for(k = m;
u.
col.
idx[k] != prow; ++k)
2020 for(--i; i >= j; --i)
2025 for(k = m;
u.
col.
idx[k] != prow; ++k)
2085 len =
u.
row.
len[rw] + beg - 1;
2095 for(i = len - 1; i >= beg; --i)
2104 l_maxabs *= threshold;
2110 for(i = len; i >= beg; --i)
2116 if(j < mkwtz &&
spxAbs(x) > l_maxabs)
2131 len =
u.
col.
len[cl] + beg - 1;
2139 for(i = len; i >= beg; --i)
2172 for(kk = m +
u.
row.
len[k] - 1; kk >= m; --kk)
2185 for(--kk; kk > m; --kk)
2194 l_maxabs *= threshold;
2231 if(pr->
idx == rw || pr->
mkwtz >= mkwtz)
2245 if(num >= candidates)
2277 int c, i, j, k, ll, m, n;
2286 for(j = m;
u.
row.
idx[j] != pcol; --j)
2307 for(j = m - 1; j >= n; --j)
2336 for(i = k;
u.
col.
idx[i] != r; --i)
2406 int i, j, k, m = -1;
2411 int plen = --(
u.
row.
len[prow]);
2412 int pend = pbeg + plen;
2441 for(i = pbeg; i < pend; ++i)
2449 for(k = m;
u.
col.
idx[k] != prow; ++k)
2474 for(++i; i < m; ++i)
2486 for(i =
u.
row.
start[prow], pend = i + plen; i < pend; ++i)
2520 for(i = 0; i <
thedim; ++i)
2548 "ERROR: no pivot element selected");
2551 pivot = pivot->
next)
2574 "ERROR: one row must be left");
2576 "ERROR: one col must be left");
2594 for(i = 0; i <
thedim; i++)
2599 for(i = 0; i < thedim; i++)
2610 assert((*idx >= 0) && (*idx < thedim));
2614 assert((k >= 0) && (k <
u.
col.
size));
2696 for(i =
thedim; i--; *l_rbeg++ = 0)
2698 *rorig++ = *rrorig++;
2699 *rperm++ = *rrperm++;
2704 l_rbeg =
l.
rbeg + 1;
2711 for(m = 0, i =
thedim; i--; l_rbeg++)
2720 l_rbeg =
l.
rbeg + 1;
2722 for(i = j = 0; i < vecs; ++i)
2727 for(; j < beg[i + 1]; j++)
2729 k = l_rbeg[*idx++]++;
2732 l_rval[k] = val[validx];
2739 assert(
l.
rbeg[0] == 0);
2751 MSG_DEBUG(std::cout <<
"CLUFactorRational::factor()\n");
2823 for(i = 0; i <
thedim; ++i)
2826 std::cout <<
"DCLUFA01 diag[" << i <<
"]: [" <<
col.
orig[
row.
perm[i]]
2827 <<
"] = " <<
diag[i] << std::endl;
2829 for(j = 0; j <
u.
row.
len[i]; ++j)
2830 std::cout <<
"DCLUFA02 u[" << i <<
"]: [" 2837 for(i = 0; i <
thedim; ++i)
2842 std::cout <<
"DCLUFA03 l[" << i <<
"]" << std::endl;
2845 std::cout <<
"DCLUFA04 l[" << k -
l.
start[j]
2846 <<
"]: [" <<
l.
idx[k]
2847 <<
"] = " <<
l.
val[k] << std::endl;
2901 int newsize = int(0.2 *
l.
val.
dim() + size);
2917 int* p_lrow =
l.
row;
2922 assert(p_len > 0 &&
"ERROR: no empty columns allowed in L vectors");
2938 #ifdef ENABLE_CONSISTENCY_CHECKS 2955 assert(ring->
idx >= 0);
2980 assert(ring->
idx >= 0);
2997 for(i = 0; i <
thedim; ++i)
3006 for(i = 0; i <
thedim; ++i)
3033 for(i = 0; i <
thedim; ++i)
3055 for(i = 0; i < thedim -
temp.
stage; ++i)
3069 #endif // CONSISTENCY_CHECKS 3077 for(
int i =
thedim - 1; i >= 0; i--)
3101 int* cidx, *clen, *cbeg;
3117 for(i =
thedim - 1; i >= 0; --i)
3120 x =
diag[r] * rhs[r];
3127 val = &cval[cbeg[c]];
3128 idx = &cidx[cbeg[c]];
3132 rhs[*idx++] -= x * (*val++);
3144 int* cidx, *clen, *cbeg;
3158 for(i =
thedim - 1; i >= 0; --i)
3162 p_work1[c] = x1 =
diag[r] * vec1[r];
3163 p_work2[c] = x2 =
diag[r] * vec2[r];
3164 vec1[r] = vec2[r] = 0;
3166 if(x1 != 0 && x2 != 0)
3168 val = &cval[cbeg[c]];
3169 idx = &cidx[cbeg[c]];
3174 vec1[*idx] -= x1 * (*val);
3175 vec2[*idx++] -= x2 * (*val++);
3180 val = &cval[cbeg[c]];
3181 idx = &cidx[cbeg[c]];
3185 vec1[*idx++] -= x1 * (*val++);
3189 val = &cval[cbeg[c]];
3190 idx = &cidx[cbeg[c]];
3194 vec2[*idx++] -= x2 * (*val++);
3204 int* cidx, *clen, *cbeg;
3205 bool notzero1, notzero2;
3221 for(i =
thedim - 1; i >= 0; --i)
3225 p_work1[c] = x1 =
diag[r] * vec1[r];
3226 p_work2[c] = x2 =
diag[r] * vec2[r];
3227 vec1[r] = vec2[r] = 0;
3228 notzero1 = x1 != 0 ? 1 : 0;
3229 notzero2 = x2 != 0 ? 1 : 0;
3231 if(notzero1 && notzero2)
3235 val = &cval[cbeg[c]];
3236 idx = &cidx[cbeg[c]];
3241 vec1[*idx] -= x1 * (*val);
3242 vec2[*idx++] -= x2 * (*val++);
3250 val = &cval[cbeg[c]];
3251 idx = &cidx[cbeg[c]];
3255 vec1[*idx++] -= x1 * (*val++);
3260 val = &cval[cbeg[c]];
3261 idx = &cidx[cbeg[c]];
3265 vec2[*idx++] -= x2 * (*val++);
3283 int* lrow, *lidx, *idx;
3293 for(i = 0; i < end; ++i)
3295 if((x = vec[lrow[i]]) != 0)
3300 MSG_DEBUG(std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl;)
3306 for(j = lbeg[i + 1]; j > k; --j)
3308 MSG_DEBUG(std::cout <<
" -> y" << *idx <<
" -= " << x <<
" * " << *val <<
3309 " = " << x * (*val) <<
" -> " << vec[*idx] - x * (*val) << std::endl;)
3310 vec[*idx++] -= x * (*val++);
3317 MSG_DEBUG(std::cout <<
"performing FT updates..." << std::endl;)
3328 for(j = lbeg[i + 1]; j > k; --j)
3329 x += vec[*idx++] * (*val++);
3333 MSG_DEBUG(std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl;)
3336 MSG_DEBUG(std::cout <<
"finished FT updates." << std::endl;)
3347 int* lrow, *lidx, *idx;
3357 for(i = 0; i < end; ++i)
3362 if(x1 != 0 && x2 != 0)
3368 for(j = lbeg[i + 1]; j > k; --j)
3370 vec1[*idx] -= x1 * (*val);
3371 vec2[*idx++] -= x2 * (*val++);
3380 for(j = lbeg[i + 1]; j > k; --j)
3381 vec1[*idx++] -= x1 * (*val++);
3389 for(j = lbeg[i + 1]; j > k; --j)
3390 vec2[*idx++] -= x2 * (*val++);
3406 for(j = lbeg[i + 1]; j > k; --j)
3408 x1 += vec1[*idx] * (*val);
3409 x2 += vec2[*idx++] * (*val++);
3412 vec1[lrow[i]] -= x1;
3414 vec2[lrow[i]] -= x2;
3425 int* lrow, *lidx, *idx;
3439 if((x = vec[lrow[i]]) != 0)
3445 for(j = lbeg[i + 1]; j > k; --j)
3446 vec[*idx++] -= x * (*val++);
3476 if(x1 != 0 && x2 != 0)
3482 for(j = lbeg[i + 1]; j > k; --j)
3484 vec1[*idx] -= x1 * (*val);
3485 vec2[*idx++] -= x2 * (*val++);
3494 for(j = lbeg[i + 1]; j > k; --j)
3495 vec1[*idx++] -= x1 * (*val++);
3503 for(j = lbeg[i + 1]; j > k; --j)
3504 vec2[*idx++] -= x2 * (*val++);
3518 for(
int i = 0; i <
thedim; i++)
3522 n += rhs[i] != 0 ? 1 : 0;
3562 for(
int i = 0; i <
thedim; i++)
3565 forest[i] = rhs1[i];
3566 n += rhs1[i] != 0 ? 1 : 0;
3602 for(
int i = 0; i <
thedim; ++i)
3631 for(
int m =
u.
row.
start[r]; m < end; m++)
3653 int* ridx, *rlen, *rbeg, *idx;
3664 for(i = 0; i <
thedim; ++i)
3671 if((x1 != 0) && (x2 != 0))
3681 for(
int m = rlen[r]; m != 0; --m)
3683 vec1[*idx] -= x1 * (*val);
3684 vec2[*idx] -= x2 * (*val++);
3696 for(
int m = rlen[r]; m != 0; --m)
3697 vec1[*idx++] -= x1 * (*val++);
3707 for(
int m = rlen[r]; m != 0; --m)
3708 vec2[*idx++] -= x2 * (*val++);
3721 int* lidx, *idx, *lrow;
3745 for(j = lbeg[i + 1]; j > k; --j)
3747 vec1[*idx] -= x1 * (*val);
3748 vec2[*idx++] -= x2 * (*val++);
3757 for(j = lbeg[i + 1]; j > k; --j)
3758 vec1[*idx++] -= x1 * (*val++);
3767 for(j = lbeg[i + 1]; j > k; --j)
3768 vec2[*idx++] -= x2 * (*val++);
3807 for(j = lbeg[i + 1]; j > k; --j)
3809 x1 += vec1[*idx] * (*val);
3810 x2 += vec2[*idx++] * (*val++);
3813 vec1[lrow[i]] -= x1;
3815 vec2[lrow[i]] -= x2;
3828 if(x1not0 && x2not0)
3831 j = rbeg[r + 1] - k;
3838 vec1[*idx] -= x1 * *val;
3839 vec2[*idx++] -= x2 * *val++;
3845 j = rbeg[r + 1] - k;
3852 vec1[*idx++] -= x1 * *val++;
3858 j = rbeg[r + 1] - k;
3865 vec2[*idx++] -= x2 * *val++;
3878 int* idx, *lidx, *lrow, *lbeg;
3889 if((x = vec[lrow[i]]) != 0)
3898 for(j = lbeg[i + 1]; j > k; --j)
3899 vec[*idx++] -= x * (*val++);
3927 for(
int j = lbeg[i + 1]; j > k; --j)
3928 x += vec[*idx++] * (*val++);
3935 for(
int i =
thedim - 1; i >= 0; --i)
3945 for(
int k =
l.
rbeg[r]; k <
l.
rbeg[r + 1]; k++)
3951 vec[j] -= x *
l.
rval[k];
3956 #endif // WITH_L_ROWS 3987 for(j = lbeg[i + 1]; j > k; --j)
3988 x += vec[*idx++] * (*val++);
4005 j = rbeg[r + 1] - k;
4012 vec[*idx++] -= x * *val++;
4029 int* lrow, *lidx, *idx;
4048 for(j = lbeg[i + 1]; j > k; --j)
4049 x += vec[*idx++] * (*val++);
4060 int* lrow, *lidx, *idx;
4080 for(j = lbeg[i + 1]; j > k; --j)
4082 x1 += vec1[*idx] * (*val);
4083 x2 += vec2[*idx++] * (*val++);
4086 vec1[lrow[i]] -= x1;
4088 vec2[lrow[i]] -= x2;
4097 int* lrow, *lidx, *idx;
4112 assert(k >= 0 && k <
l.
val.
dim());
4117 for(j = lbeg[i + 1]; j > k; --j)
4119 assert(*idx >= 0 && *idx <
thedim);
4120 x += vec[*idx++] * (*val++);
4207 Rational* rhs,
int* rhsidx,
int rhsn)
4210 int i, j, k, n, r, c;
4211 int* rorig, *corig, *cperm;
4212 int* ridx, *rlen, *rbeg, *idx;
4222 for(i = 0; i < rhsn;)
4238 assert(i >= 0 && i <
thedim);
4240 assert(c >= 0 && c <
thedim);
4247 assert(r >= 0 && r <
thedim);
4256 for(
int m = rlen[r]; m; --m)
4259 assert(j >= 0 && j <
thedim);
4290 int* rorig, *corig, *cperm;
4291 int* ridx, *rlen, *rbeg, *idx;
4301 for(i = 0; i < rhsn;)
4315 assert(i >= 0 && i <
thedim);
4317 assert(c >= 0 && c <
thedim);
4324 assert(r >= 0 && r <
thedim);
4332 for(
int m = rlen[r]; m; --m)
4335 assert(j >= 0 && j <
thedim);
4365 int* idx, *lidx, *lrow, *lbeg;
4375 assert(i >= 0 && i <
l.
val.
dim());
4377 if((x = vec[lrow[i]]) != 0)
4380 assert(k >= 0 && k <
l.
val.
dim());
4384 for(j = lbeg[i + 1]; j > k; --j)
4387 assert(m >= 0 && m <
thedim);
4418 int* idx, *lidx, *lrow, *lbeg;
4428 if((x = vec[lrow[i]]) != 0)
4430 assert(i >= 0 && i <
l.
val.
dim());
4432 assert(k >= 0 && k <
l.
val.
dim());
4436 for(j = lbeg[i + 1]; j > k; --j)
4438 assert(*idx >= 0 && *idx <
thedim);
4439 vec[*idx++] -= x * (*val++);
4466 #pragma warn "Not yet implemented, define WITH_L_ROWS" 4479 for(j = lbeg[i + 1]; j > k; --j)
4480 x += vec[*idx++] * (*val++);
4505 j = rbeg[r + 1] - k;
4511 assert(
l.
rperm[*idx] < i);
4537 for(i = 0; i < n; ++i)
4573 assert(k >= 0 && k <
l.
val.
dim());
4578 for(j = lbeg[i + 1]; j > k; --j)
4580 assert(*idx >= 0 && *idx <
thedim);
4581 x += vec[*idx++] * (*val++);
4597 j = rbeg[r + 1] - k;
4603 assert(
l.
rperm[*idx] < i);
4604 vec[*idx++] -= x * *val++;
4618 int* lrow, *lidx, *idx;
4628 for(i = 0; i < end; ++i)
4638 for(j = lbeg[i + 1]; j > k; --j)
4640 assert(*idx >= 0 && *idx <
thedim);
4641 ridx[rn] = n = *idx++;
4642 rn += (vec[n] == 0) ? 1 : 0;
4643 vec[n] -= x * (*val++);
4660 for(j = lbeg[i + 1]; j > k; --j)
4662 assert(*idx >= 0 && *idx <
thedim);
4663 x += vec[*idx++] * (*val++);
4666 ridx[rn] = j = lrow[i];
4668 rn += (vec[j] == 0) ? 1 : 0;
4678 Rational* vec2,
int* ridx2,
int* rn2ptr)
4685 int* lrow, *lidx, *idx;
4698 for(i = 0; i < end; ++i)
4712 for(j = lbeg[i + 1]; j > k; --j)
4714 assert(*idx >= 0 && *idx <
thedim);
4715 ridx[rn] = ridx2[rn2] = n = *idx++;
4718 rn += (y == 0) ? 1 : 0;
4719 rn2 += (y2 == 0) ? 1 : 0;
4721 y2 -= x2 * (*val++);
4730 for(j = lbeg[i + 1]; j > k; --j)
4732 assert(*idx >= 0 && *idx <
thedim);
4733 ridx[rn] = n = *idx++;
4735 rn += (y == 0) ? 1 : 0;
4748 for(j = lbeg[i + 1]; j > k; --j)
4750 assert(*idx >= 0 && *idx <
thedim);
4751 ridx2[rn2] = n = *idx++;
4753 rn2 += (y2 == 0) ? 1 : 0;
4754 y2 -= x2 * (*val++);
4772 for(j = lbeg[i + 1]; j > k; --j)
4774 assert(*idx >= 0 && *idx <
thedim);
4775 x += vec[*idx] * (*val);
4776 x2 += vec2[*idx++] * (*val++);
4779 ridx[rn] = ridx2[rn2] = j = lrow[i];
4781 rn += (vec[j] == 0) ? 1 : 0;
4782 rn2 += (vec2[j] == 0) ? 1 : 0;
4796 Rational* vec2,
int* ridx2,
int* rn2ptr,
4797 Rational* vec3,
int* ridx3,
int* rn3ptr)
4805 int* lrow, *lidx, *idx;
4819 for(i = 0; i < end; ++i)
4837 for(j = lbeg[i + 1]; j > k; --j)
4839 assert(*idx >= 0 && *idx <
thedim);
4840 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4844 rn += (y == 0) ? 1 : 0;
4845 rn2 += (y2 == 0) ? 1 : 0;
4846 rn3 += (y3 == 0) ? 1 : 0;
4849 y3 -= x3 * (*val++);
4861 for(j = lbeg[i + 1]; j > k; --j)
4863 assert(*idx >= 0 && *idx <
thedim);
4864 ridx[rn] = ridx2[rn2] = n = *idx++;
4867 rn += (y == 0) ? 1 : 0;
4868 rn2 += (y2 == 0) ? 1 : 0;
4870 y2 -= x2 * (*val++);
4881 for(j = lbeg[i + 1]; j > k; --j)
4883 assert(*idx >= 0 && *idx <
thedim);
4884 ridx[rn] = ridx3[rn3] = n = *idx++;
4887 rn += (y == 0) ? 1 : 0;
4888 rn3 += (y3 == 0) ? 1 : 0;
4890 y3 -= x3 * (*val++);
4900 for(j = lbeg[i + 1]; j > k; --j)
4902 assert(*idx >= 0 && *idx <
thedim);
4903 ridx[rn] = n = *idx++;
4905 rn += (y == 0) ? 1 : 0;
4921 for(j = lbeg[i + 1]; j > k; --j)
4923 assert(*idx >= 0 && *idx <
thedim);
4924 ridx2[rn2] = ridx3[rn3] = n = *idx++;
4927 rn2 += (y2 == 0) ? 1 : 0;
4928 rn3 += (y3 == 0) ? 1 : 0;
4930 y3 -= x3 * (*val++);
4940 for(j = lbeg[i + 1]; j > k; --j)
4942 assert(*idx >= 0 && *idx <
thedim);
4943 ridx2[rn2] = n = *idx++;
4945 rn2 += (y2 == 0) ? 1 : 0;
4946 y2 -= x2 * (*val++);
4959 for(j = lbeg[i + 1]; j > k; --j)
4961 assert(*idx >= 0 && *idx <
thedim);
4962 ridx3[rn3] = n = *idx++;
4964 rn3 += (y3 == 0) ? 1 : 0;
4965 y3 -= x3 * (*val++);
4983 for(j = lbeg[i + 1]; j > k; --j)
4985 assert(*idx >= 0 && *idx <
thedim);
4986 x += vec[*idx] * (*val);
4987 x2 += vec2[*idx] * (*val);
4988 x3 += vec3[*idx++] * (*val++);
4991 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
4993 rn += (vec[j] == 0) ? 1 : 0;
4994 rn2 += (vec2[j] == 0) ? 1 : 0;
4995 rn3 += (vec3[j] == 0) ? 1 : 0;
5014 int i, j, k, r, c, n;
5017 int* cidx, *clen, *cbeg;
5039 assert(i >= 0 && i <
thedim);
5041 assert(r >= 0 && r <
thedim);
5043 x =
diag[r] * rhs[r];
5049 assert(c >= 0 && c <
thedim);
5052 val = &cval[cbeg[c]];
5053 idx = &cidx[cbeg[c]];
5058 assert(*idx >= 0 && *idx <
thedim);
5060 assert(k >= 0 && k <
thedim);
5081 if(rn > i * verySparseFactor4right)
5084 for(i = *ridx; i >= 0; --i)
5087 assert(r >= 0 && r <
thedim);
5088 x =
diag[r] * rhs[r];
5094 assert(c >= 0 && c <
thedim);
5097 val = &cval[cbeg[c]];
5098 idx = &cidx[cbeg[c]];
5103 assert(*idx >= 0 && *idx <
thedim);
5104 rhs[*idx++] -= x * (*val++);
5123 int* cidx, *clen, *cbeg;
5140 if(rn > *ridx * verySparseFactor4right)
5143 for(i = *ridx; i >= 0; --i)
5145 assert(i >= 0 && i <
thedim);
5147 assert(r >= 0 && r <
thedim);
5148 x =
diag[r] * rhs[r];
5155 val = &cval[cbeg[c]];
5156 idx = &cidx[cbeg[c]];
5161 assert(*idx >= 0 && *idx <
thedim);
5162 rhs[*idx++] -= x * (*val++);
5174 assert(i >= 0 && i <
thedim);
5178 assert(r >= 0 && r <
thedim);
5180 x =
diag[r] * rhs[r];
5188 val = &cval[cbeg[c]];
5189 idx = &cidx[cbeg[c]];
5195 assert(k >= 0 && k <
thedim);
5223 int i, j, k, r, c, n;
5226 int* cidx, *clen, *cbeg;
5252 else if(*ridx2 > *ridx)
5254 else if(*ridx2 < *ridx)
5262 assert(i >= 0 && i <
thedim);
5265 assert(r >= 0 && r <
thedim);
5267 x =
diag[r] * rhs[r];
5268 x2 =
diag[r] * rhs2[r];
5278 val = &cval[cbeg[c]];
5279 idx = &cidx[cbeg[c]];
5287 assert(k >= 0 && k <
thedim);
5332 assert(k >= 0 && k <
thedim);
5357 assert(c >= 0 && c <
thedim);
5359 val = &cval[cbeg[c]];
5360 idx = &cidx[cbeg[c]];
5366 assert(k >= 0 && k <
thedim);
5371 y2 = -x2 * (*val++);
5381 y2 -= x2 * (*val++);
5388 if(rn + rn2 > i * verySparseFactor4right)
5400 assert(r >= 0 && r <
thedim);
5401 x =
diag[r] * rhs[r];
5402 x2 =
diag[r] * rhs2[r];
5409 assert(c >= 0 && c <
thedim);
5411 val = &cval[cbeg[c]];
5412 idx = &cidx[cbeg[c]];
5422 assert(*idx >= 0 && *idx <
thedim);
5423 rhs[*idx] -= x * (*val);
5424 rhs2[*idx++] -= x2 * (*val++);
5431 assert(*idx >= 0 && *idx <
thedim);
5432 rhs2[*idx++] -= x2 * (*val++);
5439 assert(c >= 0 && c <
thedim);
5442 val = &cval[cbeg[c]];
5443 idx = &cidx[cbeg[c]];
5448 assert(*idx >= 0 && *idx <
thedim);
5449 rhs[*idx++] -= x * (*val++);
5467 int* lrow, *lidx, *idx;
5480 assert(i >= 0 && i <
thedim);
5486 assert(k >= 0 && k <
l.
val.
dim());
5490 for(j = lbeg[i + 1]; j > k; --j)
5492 int m = ridx[n] = *idx++;
5493 assert(m >= 0 && m <
thedim);
5495 n += (y == 0) ? 1 : 0;
5496 y = y - x * (*val++);
5512 int* lrow, *lidx, *idx;
5525 assert(i >= 0 && i <
thedim);
5527 if((x = vec[lrow[i]]) != 0)
5530 assert(k >= 0 && k <
l.
val.
dim());
5534 for(j = lbeg[i + 1]; j > k; --j)
5536 assert(*idx >= 0 && *idx <
thedim);
5537 vec[*idx++] -= x * (*val++);
5547 Rational* forest,
int* forestNum,
int* forestIdx)
5559 int* it = forestIdx;
5563 for(i = j = 0; i < rn; ++i)
5566 assert(k >= 0 && k <
thedim);
5578 *forestNum = rn = j;
5588 for(i = j = 0; i < rn; ++i)
5591 assert(k >= 0 && k <
thedim);
5615 Rational* rhs2,
int* ridx2,
int rn2,
5616 Rational* forest,
int* forestNum,
int* forestIdx)
5628 int* it = forestIdx;
5632 for(i = j = 0; i < rn; ++i)
5635 assert(k >= 0 && k <
thedim);
5647 *forestNum = rn = j;
5657 for(i = j = 0; i < rn; ++i)
5660 assert(k >= 0 && k <
thedim);
5672 if(rn2 >
thedim * verySparseFactor4right)
5687 for(i = j = 0; i < rn2; ++i)
5690 assert(k >= 0 && k <
thedim);
5727 Rational* rhs2,
int* ridx2,
int rn2,
5729 Rational* rhs3,
int* ridx3,
int rn3,
5730 Rational* forest,
int* forestNum,
int* forestIdx)
5733 vSolveLright3(rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3);
5734 assert(rn >= 0 && rn <=
thedim);
5735 assert(rn2 >= 0 && rn2 <=
thedim);
5736 assert(rn3 >= 0 && rn3 <=
thedim);
5746 int* it = forestIdx;
5750 for(i = j = 0; i < rn; ++i)
5753 assert(k >= 0 && k <
thedim);
5765 *forestNum = rn = j;
5775 for(i = j = 0; i < rn; ++i)
5778 assert(k >= 0 && k <
thedim);
5790 if(rn2 >
thedim * verySparseFactor4right)
5802 for(i = j = 0; i < rn2; ++i)
5805 assert(k >= 0 && k <
thedim);
5819 if(rn3 >
thedim * verySparseFactor4right)
5831 for(i = j = 0; i < rn3; ++i)
5834 assert(k >= 0 && k <
thedim);
5864 Rational* rhs2,
int* ridx2,
int rn2)
5867 assert(rn2 >= 0 && rn2 <=
thedim);
5869 if(rn2 >
thedim * verySparseFactor4right)
5883 for(i = j = 0; i < rn2; ++i)
5886 assert(k >= 0 && k <
thedim);
5928 Rational* rhs2,
int* ridx2,
int rn2)
5956 Rational* rhs2,
int* ridx2,
int rn2,
5958 Rational* rhs3,
int* ridx3,
int rn3)
5989 Rational* rhs2,
int* ridx2,
int rn2)
int mkwtz
markowitz number of pivot
Rational spxAbs(const Rational &r)
Absolute.
void remaxCol(int p_col, int len)
Real lMemMult
factor of minimum Memory * number of nonzeros
VectorRational val
hold nonzero values
int * start
starting positions in val and idx
int firstUpdate
number of first update L vector
int * ridx
indices of rows of L
int solveRight2update(Rational *vec1, Rational *vec2, Rational *rhs1, Rational *rhs2, int *nonz, Rational *forest, int *forestNum, int *forestIdx)
int used
used entries of array idx
int * max
maximum available nonzeros per colunn: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
Rational * work
Working array: must always be left as 0!
int solveRight4update(Rational *vec, int *nonz, Rational *rhs, Rational *forest, int *forestNum, int *forestIdx)
Pring pivots
ring of selected pivot rows
Memory allocation routines.
int * max
maximum available nonzeros per row: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
Dring list
Double linked ringlist of vector indices in the order they appear in the column file.
void solveLleft(Rational *vec)
Real colMemMult
factor of minimum Memory * number of nonzeros
int firstUnused
number of first unused L vector
struct soplex::CLUFactorRational::U::Col col
void solveUleft(Rational *work, Rational *vec)
int size() const
Number of used indices.
VectorRational diag
Array of pivot elements.
Implementation of sparse LU factorization with Rational precision.
Exception classes for SoPlex.
int thedim
dimension of factorized matrix
#define init2DR(elem, ring)
int * row
column indices of L vectors
void forestReMaxCol(int col, int len)
int solveLleft2forest(Rational *vec1, int *, Rational *vec2)
VectorRational val
hold nonzero values: this is only initialized in the end of the factorization with DEFAULT updates...
The loaded matrix is singular.
void solveUpdateLeft(Rational *vec)
Temp temp
Temporary storage.
VectorRational val
values of L vectors
void update(int p_col, Rational *p_work, const int *p_idx, int num)
int used
used entries of arrays idx and val
void eliminateRowSingletons()
static void enQueueMin(int *heap, int *size, int elem)
void factor(const SVectorRational **vec, const Rational &threshold)
pivoting threshold
int vSolveUright2(Rational *vec, int *vidx, Rational *rhs, int *ridx, int rn, Rational *vec2, Rational *rhs2, int *ridx2, int rn2)
R & value(int n)
Reference to value of n 'th nonzero.
void solveUpdateRight2(Rational *vec1, Rational *vec2)
Wrapper for GMP type mpq_class.We wrap mpq_class so that we can replace it by a double type if GMP is...
int solveUright2eps(Rational *work1, Rational *vec1, Rational *work2, Rational *vec2, int *nonz)
int updateType
type of updates to be used.
void initFactorMatrix(const SVectorRational **vec)
void eliminatePivot(int prow, int pos)
int solveUrightEps(Rational *vec, int *nonz, Rational *rhs)
int vSolveLeft(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn)
void solveLright(Rational *vec)
int solveLleftEps(Rational *vec, int *nonz)
Rational initMaxabs
maximum abs number in initail Matrix
void solveRight(Rational *vec, Rational *rhs)
virtual void start()=0
start timer, resume accounting user, system and real time.
virtual Real stop()=0
stop timer, return accounted user time.
void spx_alloc(T &p, int n=1)
Allocate memory.
void clear()
clears the structure
int vSolveUright(Rational *vec, int *vidx, Rational *rhs, int *ridx, int rn)
static const Real verySparseFactor4right
int * start
starting positions in val and idx
int vSolveRight4update(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn, Rational *forest, int *forestNum, int *forestIdx)
int solveLeftEps(Rational *vec, Rational *rhs, int *nonz)
void solveLleftNoNZ(Rational *vec)
void solveLright2(Rational *vec1, Rational *vec2)
int pos
position of pivot column in row
void solveUright2(Rational *work1, Rational *vec1, Rational *work2, Rational *vec2)
void solveUleft2(Rational *work1, Rational *vec1, Rational *work2, Rational *vec2)
int vSolveLeft2(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn, Rational *vec2, Rational *rhs2, int *ridx2, int rn2)
int solveLeft2(Rational *vec1, int *nonz, Rational *vec2, Rational *rhs1, Rational *rhs2)
void solveUright(Rational *wrk, Rational *vec)
void setPivot(const int p_stage, const int p_col, const int p_row, const Rational &val)
int & index(int n)
Reference to index of n 'th nonzero.
Real rowMemMult
factor of minimum Memory * number of nonzeros
void solveLeft(Rational *vec, Rational *rhs)
Pring * pivot_rowNZ
lists for rows to number of nonzeros
void vSolveLeftNoNZ(Rational *vec, Rational *rhs, int *ridx, int rn)
int * s_cact
lengths of columns of active submatrix
int vSolveLright(Rational *vec, int *ridx, int rn)
void vSolveUrightNoNZ(Rational *vec, Rational *rhs, int *ridx, int rn)
Dring list
Double linked ringlist of vector indices in the order they appear in the row file.
Perm row
row permutation matrices
int * rorig
original row permutation
int size
size of array idx
void selectPivots(const Rational &threshold)
void vSolveRightNoNZ(Rational *vec2, Rational *rhs2, int *ridx2, int rn2)
Timer * factorTime
Time spent in factorizations.
int stage
stage of the structure
bool isConsistent() const
Debugging, floating point type and parameter definitions.
int * idx
array of size val.dim() storing indices of L vectors
VectorRational s_max
maximum absolute value per row (or -1)
void spx_realloc(T &p, int n)
Change amount of allocated memory.
void solveRight2(Rational *vec1, Rational *vec2, Rational *rhs1, Rational *rhs2)
void vSolveLright3(Rational *vec, int *ridx, int *rnptr, Rational *vec2, int *ridx2, int *rn2ptr, Rational *vec3, int *ridx3, int *rn3ptr)
int makeLvec(int p_len, int p_row)
Collection of dense, sparse, and semi-sparse vectors.
int dim() const
Dimension of vector.
Everything should be within this namespace.
void solveLleft2(Rational *vec1, int *, Rational *vec2)
struct soplex::CLUFactorRational::U::Row row
int * perm
perm[i] permuted index from i
void remaxRow(int p_row, int len)
int * len
used nonzeros per column vector
int * len
used nonzeros per row vectors
int * idx
hold row indices of nonzeros
int updateRow(int r, int lv, int prow, int pcol, const Rational &pval)
int vSolveUpdateRight(Rational *vec, int *ridx, int n)
int * start
starting positions in val and idx
Dring * elem
Array of ring elements.
void reDim(int newdim, const bool setZero=true)
Resets VectorBase's dimension to newdim.
int vSolveLeft3(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn, Rational *vec2, Rational *rhs2, int *ridx2, int rn2, Rational *vec3, Rational *rhs3, int *ridx3, int rn3)
int factorCount
Number of factorizations.
static int deQueueMax(int *heap, int *size)
int startSize
size of array start
void vSolveUpdateRightNoNZ(Rational *vec)
void forestUpdate(int col, Rational *work, int num, int *nonz)
Performs the Forrest-Tomlin update of the LU factorization.
void eliminateNucleus(const Rational &threshold)
int vSolveRight4update2(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn, Rational *vec2, Rational *rhs2, int *ridx2, int rn2, Rational *forest, int *forestNum, int *forestIdx)
static const Real verySparseFactor
int solveLleftForest(Rational *vec, int *)
Rational maxabs
maximum abs number in L and U
Perm col
column permutation matrices
VectorRational rval
values of rows of L
Pring * pivot_colNZ
lists for columns to number of nonzeros
static void enQueueMax(int *heap, int *size, int elem)
Exception class for status exceptions during the computationsThis class is derived from the SoPlex ex...
int * rperm
original row permutation
void solveUpdateLeft2(Rational *vec1, Rational *vec2)
Pring * pivot_row
row index handlers for Real linked list
void vSolveLright2(Rational *vec, int *ridx, int *rnptr, Rational *vec2, int *ridx2, int *rn2ptr)
void updateNoClear(int p_col, const Rational *p_work, const int *p_idx, int num)
int * rbeg
start of rows in rval and ridx
void solveUpdateRight(Rational *vec)
Dring * elem
Array of ring elements.
int * idx
array of length val.dim() to hold column indices of nonzeros in val
Pring * pivot_col
column index handlers for Real linked list
void solveLleftForestNoNZ(Rational *vec)
void eliminateColSingletons()
int idx
index of pivot row
void forestMinColMem(int size)
int * orig
orig[p] original index from p
void solveUleftNoNZ(Rational *vec, Rational *rhs, int *rhsidx, int rhsn)
int vSolveRight4update3(Rational *vec, int *idx, Rational *rhs, int *ridx, int rn, Rational *vec2, Rational *rhs2, int *ridx2, int rn2, Rational *vec3, Rational *rhs3, int *ridx3, int rn3, Rational *forest, int *forestNum, int *forestIdx)
int nzCnt
number of nonzeros in U
void spx_free(T &p)
Release memory.
void init(int p_dim)
initialization
static int deQueueMin(int *heap, int *size)
SLinSolverRational::Status stat
Status indicator.