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])
267 for(
int i = 0; i <
thedim; ++i)
285 diag[p_row] = 1.0 / val;
309 for(ring = list->
next; ring != list; ring = ring->
next)
313 if(l_rbeg[l_row] != n)
319 assert(l_rlen[l_row] <= l_rmax[l_row]);
321 l_rmax[l_row] = l_rlen[l_row];
322 j = i + l_rlen[l_row];
324 for(; i < j; ++i, ++n)
327 l_ridx[n] = l_ridx[i];
328 l_rval[n] = l_rval[i];
335 goto terminatePackRows;
340 l_rmax[l_row] = l_rlen[l_row];
367 for(ring = list->
next; ring != list; ring = ring->
next)
378 cmax[colno] = clen[colno];
391 goto terminatePackColumns;
396 cmax[colno] = clen[colno];
399 terminatePackColumns :
410 assert(
u.
row.
max[p_row] < len);
414 int delta = len -
u.
row.
max[p_row];
419 delta = len -
u.
row.
max[p_row];
429 &&
"ERROR: could not allocate memory for row file");
452 &&
"ERROR: could not allocate memory for row file");
467 for(; i < k; ++i, ++j)
496 for(ring = list->
next; ring != list; ring = ring->
next)
500 if(l_cbeg[l_col] != n)
507 l_cmax[l_col] = l_clen[l_col];
508 j = i + l_clen[l_col];
511 l_cidx[n++] = l_cidx[i];
517 goto terminatePackColumns;
522 l_cmax[l_col] = l_clen[l_col];
525 terminatePackColumns :
536 assert(
u.
col.
max[p_col] < len);
540 int delta = len -
u.
col.
max[p_col];
545 delta = len -
u.
col.
max[p_col];
555 &&
"ERROR: could not allocate memory for column file");
579 &&
"ERROR: could not allocate memory for column file");
605 assert(
u.
col.
max[p_col] < len);
609 int delta = len -
u.
col.
max[p_col];
614 delta = len -
u.
col.
max[p_col];
622 &&
"ERROR: could not allocate memory for column file");
645 &&
"ERROR: could not allocate memory for column file");
701 int i, j, k, h, m, n;
733 for(i += j - 1; i >= j; --i)
739 while(ridx[k] != p_col)
764 if(num > cmax[p_col])
775 for(j = 0; j < num; ++j)
787 assert(k - cbeg[p_col] < cmax[p_col]);
794 if(rmax[i] <= rlen[i])
801 h = rbeg[i] + (rlen[i])++;
813 nzCnt += (clen[p_col] = k - cbeg[p_col]);
830 for(i = 0; i < dim; ++i)
843 clen[p_col] = k - cbeg[p_col];
852 assert(k - cbeg[p_col] < cmax[p_col]);
859 if(rmax[i] <= rlen[i])
866 h = rbeg[i] + (rlen[i])++;
878 nzCnt += (clen[p_col] = k - cbeg[p_col]);
880 if(cbeg[p_col] + cmax[p_col] ==
u.
col.
used)
883 cmax[p_col] = clen[p_col];
899 memmove(&rorig[c], &rorig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
903 for(i = c; i <= r; ++i)
911 memmove(&corig[c], &corig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
915 for(i = c; i <= r; ++i)
927 if(i < verySparseFactor * (dim - c))
944 for(i += j - 1; i >= j; --i)
949 m = --(clen[k]) + cbeg[k];
951 for(h = m; cidx[h] != rowno; --h)
967 assert((num == 0) || (nonz != 0));
975 for(i = 0; i < num; ++i)
976 assert(p_work[corig[nonz[i]]] != 0);
986 assert(p_work[k] != 0);
990 x = p_work[k] *
diag[n];
1039 diag[rowno] = 1 / x;
1046 if(rmax[rowno] < num)
1060 for(i = 0; i < num; ++i)
1084 if(clen[j] >= cmax[j])
1091 cval[cbeg[j] + clen[j]] = x;
1093 cidx[cbeg[j] + clen[j]++] = rowno;
1097 rlen[rowno] = n - rbeg[rowno];
1103 for(i += j - 1; i >= j; --i)
1106 p_work[k] = rval[i];
1107 m = --(clen[k]) + cbeg[k];
1109 for(h = m; cidx[h] != rowno; --h)
1125 for(i = c; i < r; ++i)
1132 x = p_work[k] *
diag[n];
1146 p_work[ridx[j]] -= x * rval[j];
1169 diag[rowno] = 1 / x;
1179 for(i = r + 1; i < dim; ++i)
1180 if(p_work[corig[i]] != 0)
1197 for(i = r + 1; i < dim; ++i)
1215 if(clen[j] >= cmax[j])
1222 cval[cbeg[j] + clen[j]] = x;
1224 cidx[cbeg[j] + clen[j]++] = rowno;
1228 rlen[rowno] = n - rbeg[rowno];
1238 i = rbeg[rowno] + --(rlen[rowno]);
1239 diag[rowno] = 1 / rval[i];
1241 for(j = i = --(clen[p_col]) + cbeg[p_col]; cidx[i] != rowno; --i)
1266 assert(p_work[p_col] != 0);
1267 rezi = 1 / p_work[p_col];
1275 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1278 lval[ll] = rezi * p_work[j];
1285 lval[ll] = 1 - rezi;
1288 for(--i; i >= 0; --i)
1292 lval[ll] = x = rezi * p_work[j];
1312 assert(p_work[p_col] != 0);
1313 rezi = 1 / p_work[p_col];
1319 for(i = num - 1; (j = p_idx[i]) != p_col; --i)
1322 lval[ll] = rezi * p_work[j];
1328 lval[ll] = 1 - rezi;
1331 for(--i; i >= 0; --i)
1335 lval[ll] = x = rezi * p_work[j];
1377 Dring* rring, *lastrring;
1378 Dring* cring, *lastcring;
1388 for(
int i = 0; i <
thedim; i++)
1393 for(
int i = 0; i < thedim; i++)
1404 for(
int j = 0; j < k; ++j)
1436 lastrring->
next = rring;
1444 lastcring->
next = cring;
1448 for(
int i = 0; i <
thedim; i++)
1454 rring->
prev = lastrring;
1455 lastrring->
next = rring;
1460 cring->
prev = lastcring;
1461 lastcring->
next = cring;
1485 for(
int i = 0; i <
thedim; i++)
1495 for(
int j = 0; j < psv->
size() && nnonzeros <= 1; j++)
1497 if(psv->
value(j) != 0)
1509 else if(nnonzeros == 1)
1515 for(j = 0; psv->
value(j) == 0; j++)
1518 assert(j < psv->size());
1551 assert(nnonzeros >= 2);
1554 for(
int j = 0; j < psv->
size(); j++)
1580 assert(nnonzeros >= 2);
1601 int p_col, p_row, newrow;
1621 for(j = 0; j < len; ++j)
1630 for(k = n;
u.
col.
idx[k] != p_row; ++k)
1648 if(rperm[newrow] >= 0)
1658 for(k = n;
u.
row.
idx[k] != p_col; --k)
1693 int p_row, p_col, len, rs, lk;
1702 for(i = 0; i <
thedim; ++i)
1704 if(rperm[i] < 0 &&
u.
row.
len[i] == 1)
1743 for(j = k;
u.
row.
idx[j] != p_col; --j)
1802 for(i = 0; i <
thedim; ++i)
1887 for(; (r = idx[i]) != prow; ++i)
1894 for(j = k;
u.
row.
idx[j] != pcol; --j)
1925 assert(i < len &&
"ERROR: pivot column does not contain pivot row");
1927 for(++i; i < len; ++i)
1935 for(j = k;
u.
row.
idx[j] != pcol; --j)
1997 for(i = j; (c =
u.
row.
idx[i]) != pcol; --i)
2001 for(k = m;
u.
col.
idx[k] != prow; ++k)
2027 for(--i; i >= j; --i)
2032 for(k = m;
u.
col.
idx[k] != prow; ++k)
2092 len =
u.
row.
len[rw] + beg - 1;
2102 for(i = len - 1; i >= beg; --i)
2111 l_maxabs *= threshold;
2117 for(i = len; i >= beg; --i)
2123 if(j < mkwtz &&
spxAbs(x) > l_maxabs)
2138 len =
u.
col.
len[cl] + beg - 1;
2146 for(i = len; i >= beg; --i)
2179 for(kk = m +
u.
row.
len[k] - 1; kk >= m; --kk)
2192 for(--kk; kk > m; --kk)
2201 l_maxabs *= threshold;
2238 if(pr->
idx == rw || pr->
mkwtz >= mkwtz)
2252 if(num >= candidates)
2284 int c, i, j, k, ll, m, n;
2293 for(j = m;
u.
row.
idx[j] != pcol; --j)
2314 for(j = m - 1; j >= n; --j)
2343 for(i = k;
u.
col.
idx[i] != r; --i)
2413 int i, j, k, m = -1;
2418 int plen = --(
u.
row.
len[prow]);
2419 int pend = pbeg + plen;
2448 for(i = pbeg; i < pend; ++i)
2456 for(k = m;
u.
col.
idx[k] != prow; ++k)
2482 for(++i; i < m; ++i)
2495 for(i =
u.
row.
start[prow], pend = i + plen; i < pend; ++i)
2529 for(i = 0; i <
thedim; ++i)
2557 "ERROR: no pivot element selected");
2560 pivot = pivot->
next)
2583 "ERROR: one row must be left");
2585 "ERROR: one col must be left");
2603 for(i = 0; i <
thedim; i++)
2608 for(i = 0; i < thedim; i++)
2619 assert((*idx >= 0) && (*idx < thedim));
2623 assert((k >= 0) && (k <
u.
col.
size));
2705 for(i =
thedim; i--; *l_rbeg++ = 0)
2707 *rorig++ = *rrorig++;
2708 *rperm++ = *rrperm++;
2713 l_rbeg =
l.
rbeg + 1;
2720 for(m = 0, i =
thedim; i--; l_rbeg++)
2729 l_rbeg =
l.
rbeg + 1;
2731 for(i = j = 0; i < vecs; ++i)
2736 for(; j < beg[i + 1]; j++)
2738 k = l_rbeg[*idx++]++;
2741 l_rval[k] = val[validx];
2748 assert(
l.
rbeg[0] == 0);
2760 MSG_DEBUG(std::cout <<
"CLUFactorRational::factor()\n");
2832 for(i = 0; i <
thedim; ++i)
2835 std::cout <<
"DCLUFA01 diag[" << i <<
"]: [" <<
col.
orig[
row.
perm[i]]
2836 <<
"] = " <<
diag[i] << std::endl;
2838 for(j = 0; j <
u.
row.
len[i]; ++j)
2839 std::cout <<
"DCLUFA02 u[" << i <<
"]: [" 2846 for(i = 0; i <
thedim; ++i)
2851 std::cout <<
"DCLUFA03 l[" << i <<
"]" << std::endl;
2854 std::cout <<
"DCLUFA04 l[" << k -
l.
start[j]
2855 <<
"]: [" <<
l.
idx[k]
2856 <<
"] = " <<
l.
val[k] << std::endl;
2910 int newsize = int(0.2 *
l.
val.
dim() + size);
2926 int* p_lrow =
l.
row;
2931 assert(p_len > 0 &&
"ERROR: no empty columns allowed in L vectors");
2947 #ifdef ENABLE_CONSISTENCY_CHECKS 2964 assert(ring->
idx >= 0);
2989 assert(ring->
idx >= 0);
3006 for(i = 0; i <
thedim; ++i)
3015 for(i = 0; i <
thedim; ++i)
3042 for(i = 0; i <
thedim; ++i)
3064 for(i = 0; i < thedim -
temp.
stage; ++i)
3078 #endif // CONSISTENCY_CHECKS 3086 for(
int i =
thedim - 1; i >= 0; i--)
3110 int* cidx, *clen, *cbeg;
3126 for(i =
thedim - 1; i >= 0; --i)
3129 x =
diag[r] * rhs[r];
3136 val = &cval[cbeg[c]];
3137 idx = &cidx[cbeg[c]];
3141 rhs[*idx++] -= x * (*val++);
3153 int* cidx, *clen, *cbeg;
3167 for(i =
thedim - 1; i >= 0; --i)
3171 p_work1[c] = x1 =
diag[r] * vec1[r];
3172 p_work2[c] = x2 =
diag[r] * vec2[r];
3173 vec1[r] = vec2[r] = 0;
3175 if(x1 != 0 && x2 != 0)
3177 val = &cval[cbeg[c]];
3178 idx = &cidx[cbeg[c]];
3183 vec1[*idx] -= x1 * (*val);
3184 vec2[*idx++] -= x2 * (*val++);
3189 val = &cval[cbeg[c]];
3190 idx = &cidx[cbeg[c]];
3194 vec1[*idx++] -= x1 * (*val++);
3198 val = &cval[cbeg[c]];
3199 idx = &cidx[cbeg[c]];
3203 vec2[*idx++] -= x2 * (*val++);
3213 int* cidx, *clen, *cbeg;
3214 bool notzero1, notzero2;
3230 for(i =
thedim - 1; i >= 0; --i)
3234 p_work1[c] = x1 =
diag[r] * vec1[r];
3235 p_work2[c] = x2 =
diag[r] * vec2[r];
3236 vec1[r] = vec2[r] = 0;
3237 notzero1 = x1 != 0 ? 1 : 0;
3238 notzero2 = x2 != 0 ? 1 : 0;
3240 if(notzero1 && notzero2)
3244 val = &cval[cbeg[c]];
3245 idx = &cidx[cbeg[c]];
3250 vec1[*idx] -= x1 * (*val);
3251 vec2[*idx++] -= x2 * (*val++);
3259 val = &cval[cbeg[c]];
3260 idx = &cidx[cbeg[c]];
3264 vec1[*idx++] -= x1 * (*val++);
3269 val = &cval[cbeg[c]];
3270 idx = &cidx[cbeg[c]];
3274 vec2[*idx++] -= x2 * (*val++);
3292 int* lrow, *lidx, *idx;
3302 for(i = 0; i < end; ++i)
3304 if((x = vec[lrow[i]]) != 0)
3309 MSG_DEBUG(std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl;)
3315 for(j = lbeg[i + 1]; j > k; --j)
3317 MSG_DEBUG(std::cout <<
" -> y" << *idx <<
" -= " << x <<
" * " << *val <<
3318 " = " << x * (*val) <<
" -> " << vec[*idx] - x * (*val) << std::endl;)
3319 vec[*idx++] -= x * (*val++);
3326 MSG_DEBUG(std::cout <<
"performing FT updates..." << std::endl;)
3337 for(j = lbeg[i + 1]; j > k; --j)
3338 x += vec[*idx++] * (*val++);
3342 MSG_DEBUG(std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl;)
3345 MSG_DEBUG(std::cout <<
"finished FT updates." << std::endl;)
3356 int* lrow, *lidx, *idx;
3366 for(i = 0; i < end; ++i)
3371 if(x1 != 0 && x2 != 0)
3377 for(j = lbeg[i + 1]; j > k; --j)
3379 vec1[*idx] -= x1 * (*val);
3380 vec2[*idx++] -= x2 * (*val++);
3389 for(j = lbeg[i + 1]; j > k; --j)
3390 vec1[*idx++] -= x1 * (*val++);
3398 for(j = lbeg[i + 1]; j > k; --j)
3399 vec2[*idx++] -= x2 * (*val++);
3415 for(j = lbeg[i + 1]; j > k; --j)
3417 x1 += vec1[*idx] * (*val);
3418 x2 += vec2[*idx++] * (*val++);
3421 vec1[lrow[i]] -= x1;
3423 vec2[lrow[i]] -= x2;
3434 int* lrow, *lidx, *idx;
3448 if((x = vec[lrow[i]]) != 0)
3454 for(j = lbeg[i + 1]; j > k; --j)
3455 vec[*idx++] -= x * (*val++);
3485 if(x1 != 0 && x2 != 0)
3491 for(j = lbeg[i + 1]; j > k; --j)
3493 vec1[*idx] -= x1 * (*val);
3494 vec2[*idx++] -= x2 * (*val++);
3503 for(j = lbeg[i + 1]; j > k; --j)
3504 vec1[*idx++] -= x1 * (*val++);
3512 for(j = lbeg[i + 1]; j > k; --j)
3513 vec2[*idx++] -= x2 * (*val++);
3527 for(
int i = 0; i <
thedim; i++)
3531 n += rhs[i] != 0 ? 1 : 0;
3571 for(
int i = 0; i <
thedim; i++)
3574 forest[i] = rhs1[i];
3575 n += rhs1[i] != 0 ? 1 : 0;
3611 for(
int i = 0; i <
thedim; ++i)
3633 for(
int m =
u.
row.
start[r]; m < end; m++)
3650 int* ridx, *rlen, *rbeg, *idx;
3661 for(i = 0; i <
thedim; ++i)
3668 if((x1 != 0) && (x2 != 0))
3678 for(
int m = rlen[r]; m != 0; --m)
3680 vec1[*idx] -= x1 * (*val);
3681 vec2[*idx] -= x2 * (*val++);
3693 for(
int m = rlen[r]; m != 0; --m)
3694 vec1[*idx++] -= x1 * (*val++);
3704 for(
int m = rlen[r]; m != 0; --m)
3705 vec2[*idx++] -= x2 * (*val++);
3718 int* lidx, *idx, *lrow;
3742 for(j = lbeg[i + 1]; j > k; --j)
3744 vec1[*idx] -= x1 * (*val);
3745 vec2[*idx++] -= x2 * (*val++);
3754 for(j = lbeg[i + 1]; j > k; --j)
3755 vec1[*idx++] -= x1 * (*val++);
3764 for(j = lbeg[i + 1]; j > k; --j)
3765 vec2[*idx++] -= x2 * (*val++);
3804 for(j = lbeg[i + 1]; j > k; --j)
3806 x1 += vec1[*idx] * (*val);
3807 x2 += vec2[*idx++] * (*val++);
3810 vec1[lrow[i]] -= x1;
3812 vec2[lrow[i]] -= x2;
3825 if(x1not0 && x2not0)
3828 j = rbeg[r + 1] - k;
3835 vec1[*idx] -= x1 * *val;
3836 vec2[*idx++] -= x2 * *val++;
3842 j = rbeg[r + 1] - k;
3849 vec1[*idx++] -= x1 * *val++;
3855 j = rbeg[r + 1] - k;
3862 vec2[*idx++] -= x2 * *val++;
3875 int* idx, *lidx, *lrow, *lbeg;
3886 if((x = vec[lrow[i]]) != 0)
3895 for(j = lbeg[i + 1]; j > k; --j)
3896 vec[*idx++] -= x * (*val++);
3924 for(
int j = lbeg[i + 1]; j > k; --j)
3925 x += vec[*idx++] * (*val++);
3932 for(
int i =
thedim - 1; i >= 0; --i)
3942 for(
int k =
l.
rbeg[r]; k <
l.
rbeg[r + 1]; k++)
3948 vec[j] -= x *
l.
rval[k];
3953 #endif // WITH_L_ROWS 3984 for(j = lbeg[i + 1]; j > k; --j)
3985 x += vec[*idx++] * (*val++);
4002 j = rbeg[r + 1] - k;
4009 vec[*idx++] -= x * *val++;
4026 int* lrow, *lidx, *idx;
4045 for(j = lbeg[i + 1]; j > k; --j)
4046 x += vec[*idx++] * (*val++);
4057 int* lrow, *lidx, *idx;
4077 for(j = lbeg[i + 1]; j > k; --j)
4079 x1 += vec1[*idx] * (*val);
4080 x2 += vec2[*idx++] * (*val++);
4083 vec1[lrow[i]] -= x1;
4085 vec2[lrow[i]] -= x2;
4094 int* lrow, *lidx, *idx;
4109 assert(k >= 0 && k <
l.
val.
dim());
4114 for(j = lbeg[i + 1]; j > k; --j)
4116 assert(*idx >= 0 && *idx <
thedim);
4117 x += vec[*idx++] * (*val++);
4204 Rational* rhs,
int* rhsidx,
int rhsn)
4207 int i, j, k, n, r, c;
4208 int* rorig, *corig, *cperm;
4209 int* ridx, *rlen, *rbeg, *idx;
4219 for(i = 0; i < rhsn;)
4235 assert(i >= 0 && i <
thedim);
4237 assert(c >= 0 && c <
thedim);
4244 assert(r >= 0 && r <
thedim);
4253 for(
int m = rlen[r]; m; --m)
4256 assert(j >= 0 && j <
thedim);
4287 int* rorig, *corig, *cperm;
4288 int* ridx, *rlen, *rbeg, *idx;
4298 for(i = 0; i < rhsn;)
4312 assert(i >= 0 && i <
thedim);
4314 assert(c >= 0 && c <
thedim);
4321 assert(r >= 0 && r <
thedim);
4329 for(
int m = rlen[r]; m; --m)
4332 assert(j >= 0 && j <
thedim);
4362 int* idx, *lidx, *lrow, *lbeg;
4372 assert(i >= 0 && i <
l.
val.
dim());
4374 if((x = vec[lrow[i]]) != 0)
4377 assert(k >= 0 && k <
l.
val.
dim());
4381 for(j = lbeg[i + 1]; j > k; --j)
4384 assert(m >= 0 && m <
thedim);
4415 int* idx, *lidx, *lrow, *lbeg;
4425 if((x = vec[lrow[i]]) != 0)
4427 assert(i >= 0 && i <
l.
val.
dim());
4429 assert(k >= 0 && k <
l.
val.
dim());
4433 for(j = lbeg[i + 1]; j > k; --j)
4435 assert(*idx >= 0 && *idx <
thedim);
4436 vec[*idx++] -= x * (*val++);
4463 #pragma warn "Not yet implemented, define WITH_L_ROWS" 4476 for(j = lbeg[i + 1]; j > k; --j)
4477 x += vec[*idx++] * (*val++);
4502 j = rbeg[r + 1] - k;
4508 assert(
l.
rperm[*idx] < i);
4534 for(i = 0; i < n; ++i)
4570 assert(k >= 0 && k <
l.
val.
dim());
4575 for(j = lbeg[i + 1]; j > k; --j)
4577 assert(*idx >= 0 && *idx <
thedim);
4578 x += vec[*idx++] * (*val++);
4594 j = rbeg[r + 1] - k;
4600 assert(
l.
rperm[*idx] < i);
4601 vec[*idx++] -= x * *val++;
4615 int* lrow, *lidx, *idx;
4625 for(i = 0; i < end; ++i)
4635 for(j = lbeg[i + 1]; j > k; --j)
4637 assert(*idx >= 0 && *idx <
thedim);
4638 ridx[rn] = n = *idx++;
4639 rn += (vec[n] == 0) ? 1 : 0;
4640 vec[n] -= x * (*val++);
4657 for(j = lbeg[i + 1]; j > k; --j)
4659 assert(*idx >= 0 && *idx <
thedim);
4660 x += vec[*idx++] * (*val++);
4663 ridx[rn] = j = lrow[i];
4665 rn += (vec[j] == 0) ? 1 : 0;
4675 Rational* vec2,
int* ridx2,
int* rn2ptr)
4682 int* lrow, *lidx, *idx;
4695 for(i = 0; i < end; ++i)
4709 for(j = lbeg[i + 1]; j > k; --j)
4711 assert(*idx >= 0 && *idx <
thedim);
4712 ridx[rn] = ridx2[rn2] = n = *idx++;
4715 rn += (y == 0) ? 1 : 0;
4716 rn2 += (y2 == 0) ? 1 : 0;
4718 y2 -= x2 * (*val++);
4727 for(j = lbeg[i + 1]; j > k; --j)
4729 assert(*idx >= 0 && *idx <
thedim);
4730 ridx[rn] = n = *idx++;
4732 rn += (y == 0) ? 1 : 0;
4745 for(j = lbeg[i + 1]; j > k; --j)
4747 assert(*idx >= 0 && *idx <
thedim);
4748 ridx2[rn2] = n = *idx++;
4750 rn2 += (y2 == 0) ? 1 : 0;
4751 y2 -= x2 * (*val++);
4769 for(j = lbeg[i + 1]; j > k; --j)
4771 assert(*idx >= 0 && *idx <
thedim);
4772 x += vec[*idx] * (*val);
4773 x2 += vec2[*idx++] * (*val++);
4776 ridx[rn] = ridx2[rn2] = j = lrow[i];
4778 rn += (vec[j] == 0) ? 1 : 0;
4779 rn2 += (vec2[j] == 0) ? 1 : 0;
4793 Rational* vec2,
int* ridx2,
int* rn2ptr,
4794 Rational* vec3,
int* ridx3,
int* rn3ptr)
4802 int* lrow, *lidx, *idx;
4816 for(i = 0; i < end; ++i)
4834 for(j = lbeg[i + 1]; j > k; --j)
4836 assert(*idx >= 0 && *idx <
thedim);
4837 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4841 rn += (y == 0) ? 1 : 0;
4842 rn2 += (y2 == 0) ? 1 : 0;
4843 rn3 += (y3 == 0) ? 1 : 0;
4846 y3 -= x3 * (*val++);
4858 for(j = lbeg[i + 1]; j > k; --j)
4860 assert(*idx >= 0 && *idx <
thedim);
4861 ridx[rn] = ridx2[rn2] = n = *idx++;
4864 rn += (y == 0) ? 1 : 0;
4865 rn2 += (y2 == 0) ? 1 : 0;
4867 y2 -= x2 * (*val++);
4878 for(j = lbeg[i + 1]; j > k; --j)
4880 assert(*idx >= 0 && *idx <
thedim);
4881 ridx[rn] = ridx3[rn3] = n = *idx++;
4884 rn += (y == 0) ? 1 : 0;
4885 rn3 += (y3 == 0) ? 1 : 0;
4887 y3 -= x3 * (*val++);
4897 for(j = lbeg[i + 1]; j > k; --j)
4899 assert(*idx >= 0 && *idx <
thedim);
4900 ridx[rn] = n = *idx++;
4902 rn += (y == 0) ? 1 : 0;
4918 for(j = lbeg[i + 1]; j > k; --j)
4920 assert(*idx >= 0 && *idx <
thedim);
4921 ridx2[rn2] = ridx3[rn3] = n = *idx++;
4924 rn2 += (y2 == 0) ? 1 : 0;
4925 rn3 += (y3 == 0) ? 1 : 0;
4927 y3 -= x3 * (*val++);
4937 for(j = lbeg[i + 1]; j > k; --j)
4939 assert(*idx >= 0 && *idx <
thedim);
4940 ridx2[rn2] = n = *idx++;
4942 rn2 += (y2 == 0) ? 1 : 0;
4943 y2 -= x2 * (*val++);
4956 for(j = lbeg[i + 1]; j > k; --j)
4958 assert(*idx >= 0 && *idx <
thedim);
4959 ridx3[rn3] = n = *idx++;
4961 rn3 += (y3 == 0) ? 1 : 0;
4962 y3 -= x3 * (*val++);
4980 for(j = lbeg[i + 1]; j > k; --j)
4982 assert(*idx >= 0 && *idx <
thedim);
4983 x += vec[*idx] * (*val);
4984 x2 += vec2[*idx] * (*val);
4985 x3 += vec3[*idx++] * (*val++);
4988 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
4990 rn += (vec[j] == 0) ? 1 : 0;
4991 rn2 += (vec2[j] == 0) ? 1 : 0;
4992 rn3 += (vec3[j] == 0) ? 1 : 0;
5011 int i, j, k, r, c, n;
5014 int* cidx, *clen, *cbeg;
5036 assert(i >= 0 && i <
thedim);
5038 assert(r >= 0 && r <
thedim);
5040 x =
diag[r] * rhs[r];
5046 assert(c >= 0 && c <
thedim);
5049 val = &cval[cbeg[c]];
5050 idx = &cidx[cbeg[c]];
5055 assert(*idx >= 0 && *idx <
thedim);
5057 assert(k >= 0 && k <
thedim);
5078 if(rn > i * verySparseFactor4right)
5081 for(i = *ridx; i >= 0; --i)
5084 assert(r >= 0 && r <
thedim);
5085 x =
diag[r] * rhs[r];
5091 assert(c >= 0 && c <
thedim);
5094 val = &cval[cbeg[c]];
5095 idx = &cidx[cbeg[c]];
5100 assert(*idx >= 0 && *idx <
thedim);
5101 rhs[*idx++] -= x * (*val++);
5120 int* cidx, *clen, *cbeg;
5137 if(rn > *ridx * verySparseFactor4right)
5140 for(i = *ridx; i >= 0; --i)
5142 assert(i >= 0 && i <
thedim);
5144 assert(r >= 0 && r <
thedim);
5145 x =
diag[r] * rhs[r];
5152 val = &cval[cbeg[c]];
5153 idx = &cidx[cbeg[c]];
5158 assert(*idx >= 0 && *idx <
thedim);
5159 rhs[*idx++] -= x * (*val++);
5171 assert(i >= 0 && i <
thedim);
5175 assert(r >= 0 && r <
thedim);
5177 x =
diag[r] * rhs[r];
5185 val = &cval[cbeg[c]];
5186 idx = &cidx[cbeg[c]];
5192 assert(k >= 0 && k <
thedim);
5220 int i, j, k, r, c, n;
5223 int* cidx, *clen, *cbeg;
5249 else if(*ridx2 > *ridx)
5251 else if(*ridx2 < *ridx)
5259 assert(i >= 0 && i <
thedim);
5262 assert(r >= 0 && r <
thedim);
5264 x =
diag[r] * rhs[r];
5265 x2 =
diag[r] * rhs2[r];
5275 val = &cval[cbeg[c]];
5276 idx = &cidx[cbeg[c]];
5284 assert(k >= 0 && k <
thedim);
5329 assert(k >= 0 && k <
thedim);
5354 assert(c >= 0 && c <
thedim);
5356 val = &cval[cbeg[c]];
5357 idx = &cidx[cbeg[c]];
5363 assert(k >= 0 && k <
thedim);
5368 y2 = -x2 * (*val++);
5378 y2 -= x2 * (*val++);
5385 if(rn + rn2 > i * verySparseFactor4right)
5397 assert(r >= 0 && r <
thedim);
5398 x =
diag[r] * rhs[r];
5399 x2 =
diag[r] * rhs2[r];
5406 assert(c >= 0 && c <
thedim);
5408 val = &cval[cbeg[c]];
5409 idx = &cidx[cbeg[c]];
5419 assert(*idx >= 0 && *idx <
thedim);
5420 rhs[*idx] -= x * (*val);
5421 rhs2[*idx++] -= x2 * (*val++);
5428 assert(*idx >= 0 && *idx <
thedim);
5429 rhs2[*idx++] -= x2 * (*val++);
5436 assert(c >= 0 && c <
thedim);
5439 val = &cval[cbeg[c]];
5440 idx = &cidx[cbeg[c]];
5445 assert(*idx >= 0 && *idx <
thedim);
5446 rhs[*idx++] -= x * (*val++);
5464 int* lrow, *lidx, *idx;
5477 assert(i >= 0 && i <
thedim);
5483 assert(k >= 0 && k <
l.
val.
dim());
5487 for(j = lbeg[i + 1]; j > k; --j)
5489 int m = ridx[n] = *idx++;
5490 assert(m >= 0 && m <
thedim);
5492 n += (y == 0) ? 1 : 0;
5493 y = y - x * (*val++);
5509 int* lrow, *lidx, *idx;
5522 assert(i >= 0 && i <
thedim);
5524 if((x = vec[lrow[i]]) != 0)
5527 assert(k >= 0 && k <
l.
val.
dim());
5531 for(j = lbeg[i + 1]; j > k; --j)
5533 assert(*idx >= 0 && *idx <
thedim);
5534 vec[*idx++] -= x * (*val++);
5544 Rational* forest,
int* forestNum,
int* forestIdx)
5556 int* it = forestIdx;
5560 for(i = j = 0; i < rn; ++i)
5563 assert(k >= 0 && k <
thedim);
5575 *forestNum = rn = j;
5585 for(i = j = 0; i < rn; ++i)
5588 assert(k >= 0 && k <
thedim);
5612 Rational* rhs2,
int* ridx2,
int rn2,
5613 Rational* forest,
int* forestNum,
int* forestIdx)
5625 int* it = forestIdx;
5629 for(i = j = 0; i < rn; ++i)
5632 assert(k >= 0 && k <
thedim);
5644 *forestNum = rn = j;
5654 for(i = j = 0; i < rn; ++i)
5657 assert(k >= 0 && k <
thedim);
5669 if(rn2 >
thedim * verySparseFactor4right)
5684 for(i = j = 0; i < rn2; ++i)
5687 assert(k >= 0 && k <
thedim);
5724 Rational* rhs2,
int* ridx2,
int rn2,
5726 Rational* rhs3,
int* ridx3,
int rn3,
5727 Rational* forest,
int* forestNum,
int* forestIdx)
5730 vSolveLright3(rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3);
5731 assert(rn >= 0 && rn <=
thedim);
5732 assert(rn2 >= 0 && rn2 <=
thedim);
5733 assert(rn3 >= 0 && rn3 <=
thedim);
5743 int* it = forestIdx;
5747 for(i = j = 0; i < rn; ++i)
5750 assert(k >= 0 && k <
thedim);
5762 *forestNum = rn = j;
5772 for(i = j = 0; i < rn; ++i)
5775 assert(k >= 0 && k <
thedim);
5787 if(rn2 >
thedim * verySparseFactor4right)
5799 for(i = j = 0; i < rn2; ++i)
5802 assert(k >= 0 && k <
thedim);
5816 if(rn3 >
thedim * verySparseFactor4right)
5828 for(i = j = 0; i < rn3; ++i)
5831 assert(k >= 0 && k <
thedim);
5861 Rational* rhs2,
int* ridx2,
int rn2)
5864 assert(rn2 >= 0 && rn2 <=
thedim);
5866 if(rn2 >
thedim * verySparseFactor4right)
5880 for(i = j = 0; i < rn2; ++i)
5883 assert(k >= 0 && k <
thedim);
5925 Rational* rhs2,
int* ridx2,
int rn2)
5953 Rational* rhs2,
int* ridx2,
int rn2,
5955 Rational* rhs3,
int* ridx3,
int rn3)
5986 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)
Exception class for out of memory exceptions.This class is derived from the SoPlex exception base cla...
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.