62 for ( i = 1; i < *size; ++i )
63 for ( j = 0; j < i; ++j )
64 assert( heap[i] != heap[j] );
76 e = heap[s = --( *size )];
79 for ( j = 0, i = 1; i < s; i = 2 * j + 1 )
112 if ( i < *size && e < heap[i] )
133 if ( elem < heap[i] )
146 for ( i = 1; i < *size; ++i )
147 for ( j = 0; j < i; ++j )
148 assert( heap[i] != heap[j] );
160 e = heap[s = --( *size )];
163 for ( j = 0, i = 1; i < s; i = 2 * j + 1 )
196 if ( i < *size && e > heap[i] )
258 for (
int i = 0; i <
thedim; ++i )
276 diag[p_row] = 1.0 / val;
300 for ( ring = list->
next; ring != list; ring = ring->
next )
304 if ( l_rbeg[l_row] != n )
310 assert( l_rlen[l_row] <= l_rmax[l_row] );
312 l_rmax[l_row] = l_rlen[l_row];
313 j = i + l_rlen[l_row];
315 for ( ; i < j; ++i, ++n )
318 l_ridx[n] = l_ridx[i];
319 l_rval[n] = l_rval[i];
324 while ( ring != list );
326 goto terminatePackRows;
331 l_rmax[l_row] = l_rlen[l_row];
358 for ( ring = list->
next; ring != list; ring = ring->
next )
362 if ( cbeg[colno] != n )
369 cmax[colno] = clen[colno];
380 while ( ring != list );
382 goto terminatePackColumns;
387 cmax[colno] = clen[colno];
390 terminatePackColumns :
401 assert(
u.
row.
max[p_row] < len );
405 int delta = len -
u.
row.
max[p_row];
410 delta = len -
u.
row.
max[p_row];
420 &&
"ERROR: could not allocate memory for row file" );
443 &&
"ERROR: could not allocate memory for row file" );
458 for ( ; i < k; ++i, ++j )
487 for ( ring = list->
next; ring != list; ring = ring->
next )
491 if ( l_cbeg[l_col] != n )
498 l_cmax[l_col] = l_clen[l_col];
499 j = i + l_clen[l_col];
502 l_cidx[n++] = l_cidx[i];
506 while ( ring != list );
508 goto terminatePackColumns;
513 l_cmax[l_col] = l_clen[l_col];
516 terminatePackColumns :
527 assert(
u.
col.
max[p_col] < len );
531 int delta = len -
u.
col.
max[p_col];
536 delta = len -
u.
col.
max[p_col];
546 &&
"ERROR: could not allocate memory for column file" );
570 &&
"ERROR: could not allocate memory for column file" );
596 assert(
u.
col.
max[p_col] < len );
600 int delta = len -
u.
col.
max[p_col];
605 delta = len -
u.
col.
max[p_col];
613 &&
"ERROR: could not allocate memory for column file" );
636 &&
"ERROR: could not allocate memory for column file" );
692 int i, j, k, h, m, n;
724 for ( i += j - 1; i >= j; --i )
728 h = --( rlen[m] ) + k;
730 while ( ridx[k] != p_col )
755 if ( num > cmax[p_col] )
766 for ( j = 0; j < num; ++j )
774 if (
spxAbs( x ) > l_maxabs )
778 assert( k - cbeg[p_col] < cmax[p_col] );
785 if ( rmax[i] <= rlen[i] )
792 h = rbeg[i] + ( rlen[i] )++;
804 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
821 for ( i = 0; i < dim; ++i )
828 if (
spxAbs( x ) > l_maxabs )
834 clen[p_col] = k - cbeg[p_col];
843 assert( k - cbeg[p_col] < cmax[p_col] );
850 if ( rmax[i] <= rlen[i] )
857 h = rbeg[i] + ( rlen[i] )++;
869 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
871 if ( cbeg[p_col] + cmax[p_col] ==
u.
col.
used )
874 cmax[p_col] = clen[p_col];
890 memmove(&rorig[c], &rorig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
894 for ( i = c; i <= r; ++i )
902 memmove(&corig[c], &corig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
906 for ( i = c; i <= r; ++i )
918 if ( i < verySparseFactor * ( dim - c ) )
935 for ( i += j - 1; i >= j; --i )
940 m = --( clen[k] ) + cbeg[k];
942 for ( h = m; cidx[h] != rowno; --h )
958 assert(( num == 0 ) || ( nonz != 0 ) );
966 for ( i = 0; i < num; ++i )
967 assert( p_work[corig[nonz[i]]] != 0 );
977 assert( p_work[k] != 0 );
981 x = p_work[k] *
diag[n];
991 if (
spxAbs( x ) > l_maxabs )
1030 diag[rowno] = 1 / x;
1037 if ( rmax[rowno] < num )
1051 for ( i = 0; i < num; ++i )
1064 if (
spxAbs( x ) > l_maxabs )
1075 if ( clen[j] >= cmax[j] )
1082 cval[cbeg[j] + clen[j]] = x;
1084 cidx[cbeg[j] + clen[j]++] = rowno;
1088 rlen[rowno] = n - rbeg[rowno];
1094 for ( i += j - 1; i >= j; --i )
1097 p_work[k] = rval[i];
1098 m = --( clen[k] ) + cbeg[k];
1100 for ( h = m; cidx[h] != rowno; --h )
1116 for ( i = c; i < r; ++i )
1120 if ( p_work[k] != 0 )
1123 x = p_work[k] *
diag[n];
1129 if (
spxAbs( x ) > l_maxabs )
1136 for ( ; j < m; ++j )
1137 p_work[ridx[j]] -= x * rval[j];
1160 diag[rowno] = 1 / x;
1170 for ( i = r + 1; i < dim; ++i )
1171 if ( p_work[corig[i]] != 0 )
1174 if ( rmax[rowno] < n )
1188 for ( i = r + 1; i < dim; ++i )
1195 if (
spxAbs( x ) > l_maxabs )
1206 if ( clen[j] >= cmax[j] )
1213 cval[cbeg[j] + clen[j]] = x;
1215 cidx[cbeg[j] + clen[j]++] = rowno;
1219 rlen[rowno] = n - rbeg[rowno];
1230 i = rbeg[rowno] + --( rlen[rowno] );
1231 diag[rowno] = 1 / rval[i];
1233 for ( j = i = --( clen[p_col] ) + cbeg[p_col]; cidx[i] != rowno; --i )
1258 assert( p_work[p_col] != 0 );
1259 rezi = 1 / p_work[p_col];
1267 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1270 lval[ll] = rezi * p_work[j];
1277 lval[ll] = 1 - rezi;
1280 for ( --i; i >= 0; --i )
1284 lval[ll] = x = rezi * p_work[j];
1304 assert( p_work[p_col] != 0 );
1305 rezi = 1 / p_work[p_col];
1311 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1314 lval[ll] = rezi * p_work[j];
1320 lval[ll] = 1 - rezi;
1323 for ( --i; i >= 0; --i )
1327 lval[ll] = x = rezi * p_work[j];
1369 Dring *rring, *lastrring;
1370 Dring *cring, *lastcring;
1380 for (
int i = 0; i <
thedim; i++ )
1385 for (
int i = 0; i <
thedim; i++ )
1396 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 )
1495 if ( nnonzeros == 0 )
1503 if ( nnonzeros == 1 )
1509 for ( j = 0; psv->
value(j) == 0; j++ )
1512 assert( j < psv->size() );
1522 x = psv->
value( j );
1545 assert( nnonzeros >= 2 );
1548 for (
int j = 0; j < psv->
size(); j++ )
1550 x = psv->
value( j );
1555 k = psv->
index( j );
1574 assert( nnonzeros >= 2 );
1595 int p_col, p_row, newrow;
1611 assert( p_row >= 0 );
1615 for ( j = 0; j < len; ++j )
1624 for ( k = n;
u.
col.
idx[k] != p_row; ++k )
1642 if ( rperm[newrow] >= 0 )
1652 for ( k = n;
u.
row.
idx[k] != p_col; --k )
1688 int p_row, p_col, len, rs, lk;
1697 for ( i = 0; i <
thedim; ++i )
1699 if ( rperm[i] < 0 &&
u.
row.
len[i] == 1 )
1715 setPivot( rs, p_col, p_row, pval );
1725 i = (
u.
col.
len[p_col] -= i );
1727 for ( ; i < len; ++i )
1738 for ( j = k;
u.
row.
idx[j] != p_col; --j )
1798 for ( i = 0; i <
thedim; ++i )
1883 for ( ; ( r = idx[i] ) != prow; ++i )
1890 for ( j = k;
u.
row.
idx[j] != pcol; --j )
1921 assert( i < len &&
"ERROR: pivot column does not contain pivot row" );
1923 for ( ++i; i < len; ++i )
1931 for ( j = k;
u.
row.
idx[j] != pcol; --j )
1993 for ( i = j; ( c =
u.
row.
idx[i] ) != pcol; --i )
1997 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2023 for ( --i; i >= j; --i )
2028 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2073 if ( candidates > 4 )
2088 len =
u.
row.
len[rw] + beg - 1;
2098 for ( i = len - 1; i >= beg; --i )
2107 l_maxabs *= threshold;
2113 for ( i = len; i >= beg; --i )
2119 if ( j < mkwtz &&
spxAbs( x ) > l_maxabs )
2135 len =
u.
col.
len[cl] + beg - 1;
2143 for ( i = len; i >= beg; --i )
2176 for ( kk = m +
u.
row.
len[k] - 1; kk >= m; --kk )
2189 for ( --kk; kk > m; --kk )
2198 l_maxabs *= threshold;
2200 if (
spxAbs( x ) > l_maxabs )
2206 if ( j <= count + 1 )
2235 if ( pr->
idx == rw || pr->
mkwtz >= mkwtz )
2241 if ( pr->
idx != rw )
2249 if ( num >= candidates )
2281 int c, i, j, k, ll, m, n;
2290 for ( j = m;
u.
row.
idx[j] != pcol; --j )
2311 for ( j = m - 1; j >= n; --j )
2340 for ( i = k;
u.
col.
idx[i] != r; --i )
2353 if ( ll + fill >
u.
row.
max[r] )
2410 int i, j, k, m = -1;
2415 int plen = --(
u.
row.
len[prow] );
2416 int pend = pbeg + plen;
2445 for ( i = pbeg; i < pend; ++i )
2453 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2478 for ( ++i; i < m; ++i )
2490 for ( i =
u.
row.
start[prow], pend = i + plen; i < pend; ++i )
2524 for ( i = 0; i <
thedim; ++i )
2553 "ERROR: no pivot element selected" );
2556 pivot = pivot->
next )
2579 "ERROR: one row must be left" );
2581 "ERROR: one col must be left" );
2599 for ( i = 0; i <
thedim; i++ )
2604 for ( i = 0; i <
thedim; i++ )
2615 assert(( *idx >= 0 ) && ( *idx < thedim ) );
2619 assert(( k >= 0 ) && ( k <
u.
col.
size ) );
2701 for ( i =
thedim; i--; *l_rbeg++ = 0 )
2703 *rorig++ = *rrorig++;
2704 *rperm++ = *rrperm++;
2709 l_rbeg =
l.
rbeg + 1;
2711 for ( i = mem; i--; )
2716 for ( m = 0, i =
thedim; i--; l_rbeg++ )
2725 l_rbeg =
l.
rbeg + 1;
2727 for ( i = j = 0; i < vecs; ++i )
2732 for ( ; j < beg[i + 1]; j++ )
2734 k = l_rbeg[*idx++]++;
2737 l_rval[k] = val[validx];
2744 assert(
l.
rbeg[0] == 0 );
2756 MSG_DEBUG( std::cout <<
"CLUFactorRational::factor()\n" );
2828 for ( i = 0; i <
thedim; ++i )
2831 std::cout <<
"DCLUFA01 diag[" << i <<
"]: [" <<
col.
orig[
row.
perm[i]]
2832 <<
"] = " <<
diag[i] << std::endl;
2834 for ( j = 0; j <
u.
row.
len[i]; ++j )
2835 std::cout <<
"DCLUFA02 u[" << i <<
"]: [" 2842 for ( i = 0; i <
thedim; ++i )
2847 std::cout <<
"DCLUFA03 l[" << i <<
"]" << std::endl;
2850 std::cout <<
"DCLUFA04 l[" << k -
l.
start[j]
2851 <<
"]: [" <<
l.
idx[k]
2852 <<
"] = " <<
l.
val[k] << std::endl;
2905 int newsize = int( 0.2 *
l.
val.
dim() + size );
2921 int* p_lrow =
l.
row;
2926 assert( p_len > 0 &&
"ERROR: no empty columns allowed in L vectors" );
2942 #ifdef ENABLE_CONSISTENCY_CHECKS 2959 assert( ring->
idx >= 0 );
2961 assert( ring->
prev->
next == ring );
2984 assert( ring->
idx >= 0 );
2986 assert( ring->
prev->
next == ring );
3001 for ( i = 0; i <
thedim; ++i )
3010 for ( i = 0; i <
thedim; ++i )
3037 for ( i = 0; i <
thedim; ++i )
3059 for ( i = 0; i < thedim -
temp.
stage; ++i )
3073 #endif // CONSISTENCY_CHECKS 3081 for (
int i =
thedim - 1; i >= 0; i-- )
3105 int *cidx, *clen, *cbeg;
3121 for ( i =
thedim - 1; i >= 0; --i )
3124 x =
diag[r] * rhs[r];
3131 val = &cval[cbeg[c]];
3132 idx = &cidx[cbeg[c]];
3136 rhs[*idx++] -= x * ( *val++ );
3147 int *cidx, *clen, *cbeg;
3161 for ( i =
thedim - 1; i >= 0; --i )
3165 p_work1[c] = x1 =
diag[r] * vec1[r];
3166 p_work2[c] = x2 =
diag[r] * vec2[r];
3167 vec1[r] = vec2[r] = 0;
3169 if ( x1 != 0 && x2 != 0 )
3171 val = &cval[cbeg[c]];
3172 idx = &cidx[cbeg[c]];
3177 vec1[*idx] -= x1 * ( *val );
3178 vec2[*idx++] -= x2 * ( *val++ );
3184 val = &cval[cbeg[c]];
3185 idx = &cidx[cbeg[c]];
3189 vec1[*idx++] -= x1 * ( *val++ );
3194 val = &cval[cbeg[c]];
3195 idx = &cidx[cbeg[c]];
3199 vec2[*idx++] -= x2 * ( *val++ );
3208 int *cidx, *clen, *cbeg;
3209 bool notzero1, notzero2;
3225 for ( i =
thedim - 1; i >= 0; --i )
3229 p_work1[c] = x1 =
diag[r] * vec1[r];
3230 p_work2[c] = x2 =
diag[r] * vec2[r];
3231 vec1[r] = vec2[r] = 0;
3232 notzero1 = x1 != 0 ? 1 : 0;
3233 notzero2 = x2 != 0 ? 1 : 0;
3235 if ( notzero1 && notzero2 )
3239 val = &cval[cbeg[c]];
3240 idx = &cidx[cbeg[c]];
3245 vec1[*idx] -= x1 * ( *val );
3246 vec2[*idx++] -= x2 * ( *val++ );
3255 val = &cval[cbeg[c]];
3256 idx = &cidx[cbeg[c]];
3260 vec1[*idx++] -= x1 * ( *val++ );
3266 val = &cval[cbeg[c]];
3267 idx = &cidx[cbeg[c]];
3271 vec2[*idx++] -= x2 * ( *val++ );
3289 int *lrow, *lidx, *idx;
3299 for ( i = 0; i < end; ++i )
3301 if (( x = vec[lrow[i]] ) != 0 )
3306 MSG_DEBUG( std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl; )
3312 for ( j = lbeg[i + 1]; j > k; --j )
3314 MSG_DEBUG( std::cout <<
" -> y" << *idx <<
" -= " << x <<
" * " << *val <<
" = " << x * ( *val ) <<
" -> " << vec[*idx] - x * ( *val ) << std::endl; )
3315 vec[*idx++] -= x * ( *val++ );
3322 MSG_DEBUG( std::cout <<
"performing FT updates..." << std::endl; )
3326 for ( ; i < end; ++i )
3333 for ( j = lbeg[i + 1]; j > k; --j )
3334 x += vec[*idx++] * ( *val++ );
3338 MSG_DEBUG( std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl; )
3341 MSG_DEBUG( std::cout <<
"finished FT updates." << std::endl; )
3352 int *lrow, *lidx, *idx;
3362 for ( i = 0; i < end; ++i )
3367 if ( x1 != 0 && x2 != 0 )
3373 for ( j = lbeg[i + 1]; j > k; --j )
3375 vec1[*idx] -= x1 * ( *val );
3376 vec2[*idx++] -= x2 * ( *val++ );
3386 for ( j = lbeg[i + 1]; j > k; --j )
3387 vec1[*idx++] -= x1 * ( *val++ );
3396 for ( j = lbeg[i + 1]; j > k; --j )
3397 vec2[*idx++] -= x2 * ( *val++ );
3405 for ( ; i < end; ++i )
3413 for ( j = lbeg[i + 1]; j > k; --j )
3415 x1 += vec1[*idx] * ( *val );
3416 x2 += vec2[*idx++] * ( *val++ );
3419 vec1[lrow[i]] -= x1;
3421 vec2[lrow[i]] -= x2;
3432 int *lrow, *lidx, *idx;
3446 if (( x = vec[lrow[i]] ) != 0 )
3452 for ( j = lbeg[i + 1]; j > k; --j )
3453 vec[*idx++] -= x * ( *val++ );
3483 if ( x1 != 0 && x2 != 0 )
3489 for ( j = lbeg[i + 1]; j > k; --j )
3491 vec1[*idx] -= x1 * ( *val );
3492 vec2[*idx++] -= x2 * ( *val++ );
3502 for ( j = lbeg[i + 1]; j > k; --j )
3503 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 )
3640 for (
int m =
u.
row.
start[r]; m < end; m++ )
3661 int *ridx, *rlen, *rbeg, *idx;
3672 for ( i = 0; i <
thedim; ++i )
3679 if (( x1 != 0 ) && ( x2 != 0 ) )
3689 for (
int m = rlen[r]; m != 0; --m )
3691 vec1[*idx] -= x1 * ( *val );
3692 vec2[*idx] -= x2 * ( *val++ );
3705 for (
int m = rlen[r]; m != 0; --m )
3706 vec1[*idx++] -= x1 * ( *val++ );
3717 for (
int m = rlen[r]; m != 0; --m )
3718 vec2[*idx++] -= x2 * ( *val++ );
3731 int *lidx, *idx, *lrow;
3755 for ( j = lbeg[i + 1]; j > k; --j )
3757 vec1[*idx] -= x1 * ( *val );
3758 vec2[*idx++] -= x2 * ( *val++ );
3767 for ( j = lbeg[i + 1]; j > k; --j )
3768 vec1[*idx++] -= x1 * ( *val++ );
3778 for ( j = lbeg[i + 1]; j > k; --j )
3779 vec2[*idx++] -= x2 * ( *val++ );
3810 for ( ; i >= 0; --i )
3818 for ( j = lbeg[i + 1]; j > k; --j )
3820 x1 += vec1[*idx] * ( *val );
3821 x2 += vec2[*idx++] * ( *val++ );
3824 vec1[lrow[i]] -= x1;
3826 vec2[lrow[i]] -= x2;
3835 x1not0 = ( x1 != 0 );
3836 x2not0 = ( x2 != 0 );
3838 if ( x1not0 && x2not0 )
3841 j = rbeg[r + 1] - k;
3848 vec1[*idx] -= x1 * *val;
3849 vec2[*idx++] -= x2 * *val++;
3856 j = rbeg[r + 1] - k;
3863 vec1[*idx++] -= x1 * *val++;
3870 j = rbeg[r + 1] - k;
3877 vec2[*idx++] -= x2 * *val++;
3890 int *idx, *lidx, *lrow, *lbeg;
3901 if (( x = vec[lrow[i]] ) != 0 )
3910 for ( j = lbeg[i + 1]; j > k; --j )
3911 vec[*idx++] -= x * ( *val++ );
3939 for (
int j = lbeg[i + 1]; j > k; --j )
3940 x += vec[*idx++] * ( *val++ );
3946 for (
int i =
thedim - 1; i >= 0; --i )
3956 for (
int k =
l.
rbeg[r]; k <
l.
rbeg[r + 1]; k++ )
3960 assert(
l.
rperm[j] < i );
3962 vec[j] -= x *
l.
rval[k];
3967 #endif // WITH_L_ROWS 3998 for ( j = lbeg[i + 1]; j > k; --j )
3999 x += vec[*idx++] * ( *val++ );
4015 j = rbeg[r + 1] - k;
4022 vec[*idx++] -= x * *val++;
4039 int *lrow, *lidx, *idx;
4058 for ( j = lbeg[i + 1]; j > k; --j )
4059 x += vec[*idx++] * ( *val++ );
4070 int *lrow, *lidx, *idx;
4090 for ( j = lbeg[i + 1]; j > k; --j )
4092 x1 += vec1[*idx] * ( *val );
4093 x2 += vec2[*idx++] * ( *val++ );
4096 vec1[lrow[i]] -= x1;
4098 vec2[lrow[i]] -= x2;
4107 int *lrow, *lidx, *idx;
4122 assert( k >= 0 && k <
l.
val.
dim() );
4127 for ( j = lbeg[i + 1]; j > k; --j )
4129 assert( *idx >= 0 && *idx <
thedim );
4130 x += vec[*idx++] * ( *val++ );
4217 Rational* rhs,
int* rhsidx,
int rhsn )
4220 int i, j, k, n, r, c;
4221 int *rorig, *corig, *cperm;
4222 int *ridx, *rlen, *rbeg, *idx;
4232 for ( i = 0; i < rhsn; )
4248 assert( i >= 0 && i <
thedim );
4250 assert( c >= 0 && c <
thedim );
4257 assert( r >= 0 && r <
thedim );
4266 for (
int m = rlen[r]; m; --m )
4269 assert( j >= 0 && j <
thedim );
4274 y = -x * ( *val++ );
4284 y -= x * ( *val++ );
4300 int *rorig, *corig, *cperm;
4301 int *ridx, *rlen, *rbeg, *idx;
4311 for ( i = 0; i < rhsn; )
4325 assert( i >= 0 && i <
thedim );
4327 assert( c >= 0 && c <
thedim );
4334 assert( r >= 0 && r <
thedim );
4342 for (
int m = rlen[r]; m; --m )
4345 assert( j >= 0 && j <
thedim );
4350 y = -x * ( *val++ );
4360 y -= x * ( *val++ );
4375 int *idx, *lidx, *lrow, *lbeg;
4385 assert( i >= 0 && i <
l.
val.
dim() );
4387 if (( x = vec[lrow[i]] ) != 0 )
4390 assert( k >= 0 && k <
l.
val.
dim() );
4394 for ( j = lbeg[i + 1]; j > k; --j )
4397 assert( m >= 0 && m <
thedim );
4402 y = -x * ( *val++ );
4412 y -= x * ( *val++ );
4428 int *idx, *lidx, *lrow, *lbeg;
4438 if (( x = vec[lrow[i]] ) != 0 )
4440 assert( i >= 0 && i <
l.
val.
dim() );
4442 assert( k >= 0 && k <
l.
val.
dim() );
4446 for ( j = lbeg[i + 1]; j > k; --j )
4448 assert( *idx >= 0 && *idx <
thedim );
4449 vec[*idx++] -= x * ( *val++ );
4476 #pragma warn "Not yet implemented, define WITH_L_ROWS" 4482 for ( ; i >= 0; --i )
4489 for ( j = lbeg[i + 1]; j > k; --j )
4490 x += vec[*idx++] * ( *val++ );
4498 for ( i = 0; i < rn; )
4514 j = rbeg[r + 1] - k;
4520 assert(
l.
rperm[*idx] < i );
4546 for ( i = 0; i < n; ++i )
4579 for ( ; i >= 0; --i )
4582 assert( k >= 0 && k <
l.
val.
dim() );
4587 for ( j = lbeg[i + 1]; j > k; --j )
4589 assert( *idx >= 0 && *idx <
thedim );
4590 x += vec[*idx++] * ( *val++ );
4605 j = rbeg[r + 1] - k;
4611 assert(
l.
rperm[*idx] < i );
4612 vec[*idx++] -= x * *val++;
4626 int *lrow, *lidx, *idx;
4636 for ( i = 0; i < end; ++i )
4646 for ( j = lbeg[i + 1]; j > k; --j )
4648 assert( *idx >= 0 && *idx <
thedim );
4649 ridx[rn] = n = *idx++;
4650 rn += ( vec[n] == 0 ) ? 1 : 0;
4651 vec[n] -= x * ( *val++ );
4661 for ( ; i < end; ++i )
4668 for ( j = lbeg[i + 1]; j > k; --j )
4670 assert( *idx >= 0 && *idx <
thedim );
4671 x += vec[*idx++] * ( *val++ );
4674 ridx[rn] = j = lrow[i];
4676 rn += ( vec[j] == 0 ) ? 1 : 0;
4686 Rational* vec2,
int* ridx2,
int* rn2ptr )
4693 int *lrow, *lidx, *idx;
4706 for ( i = 0; i < end; ++i )
4720 for ( j = lbeg[i + 1]; j > k; --j )
4722 assert( *idx >= 0 && *idx <
thedim );
4723 ridx[rn] = ridx2[rn2] = n = *idx++;
4726 rn += ( y == 0 ) ? 1 : 0;
4727 rn2 += ( y2 == 0 ) ? 1 : 0;
4729 y2 -= x2 * ( *val++ );
4738 for ( j = lbeg[i + 1]; j > k; --j )
4740 assert( *idx >= 0 && *idx <
thedim );
4741 ridx[rn] = n = *idx++;
4743 rn += ( y == 0 ) ? 1 : 0;
4744 y -= x * ( *val++ );
4757 for ( j = lbeg[i + 1]; j > k; --j )
4759 assert( *idx >= 0 && *idx <
thedim );
4760 ridx2[rn2] = n = *idx++;
4762 rn2 += ( y2 == 0 ) ? 1 : 0;
4763 y2 -= x2 * ( *val++ );
4774 for ( ; i < end; ++i )
4781 for ( j = lbeg[i + 1]; j > k; --j )
4783 assert( *idx >= 0 && *idx <
thedim );
4784 x += vec[*idx] * ( *val );
4785 x2 += vec2[*idx++] * ( *val++ );
4788 ridx[rn] = ridx2[rn2] = j = lrow[i];
4790 rn += ( vec[j] == 0 ) ? 1 : 0;
4791 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
4805 Rational* vec2,
int* ridx2,
int* rn2ptr,
4806 Rational* vec3,
int* ridx3,
int* rn3ptr)
4814 int *lrow, *lidx, *idx;
4828 for ( i = 0; i < end; ++i )
4846 for ( j = lbeg[i + 1]; j > k; --j )
4848 assert( *idx >= 0 && *idx <
thedim );
4849 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4853 rn += ( y == 0 ) ? 1 : 0;
4854 rn2 += ( y2 == 0 ) ? 1 : 0;
4855 rn3 += ( y3 == 0 ) ? 1 : 0;
4857 y2 -= x2 * ( *val );
4858 y3 -= x3 * ( *val++ );
4870 for ( j = lbeg[i + 1]; j > k; --j )
4872 assert( *idx >= 0 && *idx <
thedim );
4873 ridx[rn] = ridx2[rn2] = n = *idx++;
4876 rn += ( y == 0 ) ? 1 : 0;
4877 rn2 += ( y2 == 0 ) ? 1 : 0;
4879 y2 -= x2 * ( *val++ );
4891 for ( j = lbeg[i + 1]; j > k; --j )
4893 assert( *idx >= 0 && *idx <
thedim );
4894 ridx[rn] = ridx3[rn3] = n = *idx++;
4897 rn += ( y == 0 ) ? 1 : 0;
4898 rn3 += ( y3 == 0 ) ? 1 : 0;
4900 y3 -= x3 * ( *val++ );
4910 for ( j = lbeg[i + 1]; j > k; --j )
4912 assert( *idx >= 0 && *idx <
thedim );
4913 ridx[rn] = n = *idx++;
4915 rn += ( y == 0 ) ? 1 : 0;
4916 y -= x * ( *val++ );
4932 for ( j = lbeg[i + 1]; j > k; --j )
4934 assert( *idx >= 0 && *idx <
thedim );
4935 ridx2[rn2] = ridx3[rn3] = n = *idx++;
4938 rn2 += ( y2 == 0 ) ? 1 : 0;
4939 rn3 += ( y3 == 0 ) ? 1 : 0;
4940 y2 -= x2 * ( *val );
4941 y3 -= x3 * ( *val++ );
4951 for ( j = lbeg[i + 1]; j > k; --j )
4953 assert( *idx >= 0 && *idx <
thedim );
4954 ridx2[rn2] = n = *idx++;
4956 rn2 += ( y2 == 0 ) ? 1 : 0;
4957 y2 -= x2 * ( *val++ );
4971 for ( j = lbeg[i + 1]; j > k; --j )
4973 assert( *idx >= 0 && *idx <
thedim );
4974 ridx3[rn3] = n = *idx++;
4976 rn3 += ( y3 == 0 ) ? 1 : 0;
4977 y3 -= x3 * ( *val++ );
4988 for ( ; i < end; ++i )
4995 for ( j = lbeg[i + 1]; j > k; --j )
4997 assert( *idx >= 0 && *idx <
thedim );
4998 x += vec[*idx] * ( *val );
4999 x2 += vec2[*idx] * ( *val );
5000 x3 += vec3[*idx++] * ( *val++ );
5003 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
5005 rn += ( vec[j] == 0 ) ? 1 : 0;
5006 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
5007 rn3 += ( vec3[j] == 0 ) ? 1 : 0;
5026 int i, j, k, r, c, n;
5029 int *cidx, *clen, *cbeg;
5051 assert( i >= 0 && i <
thedim );
5053 assert( r >= 0 && r <
thedim );
5055 x =
diag[r] * rhs[r];
5061 assert( c >= 0 && c <
thedim );
5064 val = &cval[cbeg[c]];
5065 idx = &cidx[cbeg[c]];
5070 assert( *idx >= 0 && *idx <
thedim );
5072 assert( k >= 0 && k <
thedim );
5077 y = -x * ( *val++ );
5087 y -= x * ( *val++ );
5093 if ( rn > i*verySparseFactor4right )
5095 for ( i = *ridx; i >= 0; --i )
5098 assert( r >= 0 && r <
thedim );
5099 x =
diag[r] * rhs[r];
5105 assert( c >= 0 && c <
thedim );
5108 val = &cval[cbeg[c]];
5109 idx = &cidx[cbeg[c]];
5114 assert( *idx >= 0 && *idx <
thedim );
5115 rhs[*idx++] -= x * ( *val++ );
5134 int *cidx, *clen, *cbeg;
5151 if ( rn > *ridx * verySparseFactor4right )
5153 for ( i = *ridx; i >= 0; --i )
5155 assert( i >= 0 && i <
thedim );
5157 assert( r >= 0 && r <
thedim );
5158 x =
diag[r] * rhs[r];
5165 val = &cval[cbeg[c]];
5166 idx = &cidx[cbeg[c]];
5171 assert( *idx >= 0 && *idx <
thedim );
5172 rhs[*idx++] -= x * ( *val++ );
5184 assert( i >= 0 && i <
thedim );
5188 assert( r >= 0 && r <
thedim );
5190 x =
diag[r] * rhs[r];
5198 val = &cval[cbeg[c]];
5199 idx = &cidx[cbeg[c]];
5205 assert( k >= 0 && k <
thedim );
5210 y = -x * ( *val++ );
5220 y -= x * ( *val++ );
5233 int i, j, k, r, c, n;
5236 int *cidx, *clen, *cbeg;
5254 while ( rn + rn2 > 0 )
5264 if ( *ridx2 > *ridx )
5267 if ( *ridx2 < *ridx )
5275 assert( i >= 0 && i <
thedim );
5278 assert( r >= 0 && r <
thedim );
5280 x =
diag[r] * rhs[r];
5281 x2 =
diag[r] * rhs2[r];
5291 val = &cval[cbeg[c]];
5292 idx = &cidx[cbeg[c]];
5300 assert( k >= 0 && k <
thedim );
5305 y2 = -x2 * ( *val );
5315 y2 -= x2 * ( *val );
5324 y = -x * ( *val++ );
5334 y -= x * ( *val++ );
5345 assert( k >= 0 && k <
thedim );
5350 y = -x * ( *val++ );
5360 y -= x * ( *val++ );
5371 assert( c >= 0 && c <
thedim );
5373 val = &cval[cbeg[c]];
5374 idx = &cidx[cbeg[c]];
5380 assert( k >= 0 && k <
thedim );
5385 y2 = -x2 * ( *val++ );
5395 y2 -= x2 * ( *val++ );
5402 if ( rn + rn2 > i*verySparseFactor4right )
5404 if ( *ridx > *ridx2 )
5409 for ( ; i >= 0; --i )
5413 assert( r >= 0 && r <
thedim );
5414 x =
diag[r] * rhs[r];
5415 x2 =
diag[r] * rhs2[r];
5422 assert( c >= 0 && c <
thedim );
5424 val = &cval[cbeg[c]];
5425 idx = &cidx[cbeg[c]];
5435 assert( *idx >= 0 && *idx <
thedim );
5436 rhs[*idx] -= x * ( *val );
5437 rhs2[*idx++] -= x2 * ( *val++ );
5444 assert( *idx >= 0 && *idx <
thedim );
5445 rhs2[*idx++] -= x2 * ( *val++ );
5453 assert( c >= 0 && c <
thedim );
5456 val = &cval[cbeg[c]];
5457 idx = &cidx[cbeg[c]];
5462 assert( *idx >= 0 && *idx <
thedim );
5463 rhs[*idx++] -= x * ( *val++ );
5481 int *lrow, *lidx, *idx;
5494 assert( i >= 0 && i <
thedim );
5500 assert( k >= 0 && k <
l.
val.
dim() );
5504 for ( j = lbeg[i + 1]; j > k; --j )
5506 int m = ridx[n] = *idx++;
5507 assert( m >= 0 && m <
thedim );
5509 n += ( y == 0 ) ? 1 : 0;
5510 y = y - x * ( *val++ );
5526 int *lrow, *lidx, *idx;
5539 assert( i >= 0 && i <
thedim );
5541 if (( x = vec[lrow[i]] ) != 0 )
5544 assert( k >= 0 && k <
l.
val.
dim() );
5548 for ( j = lbeg[i + 1]; j > k; --j )
5550 assert( *idx >= 0 && *idx <
thedim );
5551 vec[*idx++] -= x * ( *val++ );
5560 Rational* forest,
int* forestNum,
int* forestIdx )
5572 int* it = forestIdx;
5576 for ( i = j = 0; i < rn; ++i )
5579 assert( k >= 0 && k <
thedim );
5591 *forestNum = rn = j;
5601 for ( i = j = 0; i < rn; ++i )
5604 assert( k >= 0 && k <
thedim );
5627 Rational* rhs2,
int* ridx2,
int rn2,
5628 Rational* forest,
int* forestNum,
int* forestIdx )
5640 int* it = forestIdx;
5644 for ( i = j = 0; i < rn; ++i )
5647 assert( k >= 0 && k <
thedim );
5659 *forestNum = rn = j;
5669 for ( i = j = 0; i < rn; ++i )
5672 assert( k >= 0 && k <
thedim );
5684 if ( rn2 >
thedim*verySparseFactor4right )
5699 for ( i = j = 0; i < rn2; ++i )
5702 assert( k >= 0 && k <
thedim );
5738 Rational* rhs2,
int* ridx2,
int rn2,
5740 Rational* rhs3,
int* ridx3,
int rn3,
5741 Rational* forest,
int* forestNum,
int* forestIdx )
5744 vSolveLright3( rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3 );
5745 assert( rn >= 0 && rn <=
thedim );
5746 assert( rn2 >= 0 && rn2 <=
thedim );
5747 assert( rn3 >= 0 && rn3 <=
thedim );
5757 int* it = forestIdx;
5761 for ( i = j = 0; i < rn; ++i )
5764 assert( k >= 0 && k <
thedim );
5776 *forestNum = rn = j;
5786 for ( i = j = 0; i < rn; ++i )
5789 assert( k >= 0 && k <
thedim );
5801 if ( rn2 >
thedim*verySparseFactor4right )
5813 for ( i = j = 0; i < rn2; ++i )
5816 assert( k >= 0 && k <
thedim );
5830 if ( rn3 >
thedim*verySparseFactor4right )
5842 for ( i = j = 0; i < rn3; ++i )
5845 assert( k >= 0 && k <
thedim );
5875 Rational* rhs2,
int* ridx2,
int rn2 )
5878 assert( rn2 >= 0 && rn2 <=
thedim );
5880 if ( rn2 >
thedim*verySparseFactor4right )
5894 for ( i = j = 0; i < rn2; ++i )
5897 assert( k >= 0 && k <
thedim );
5939 Rational* rhs2,
int* ridx2,
int rn2 )
5967 Rational* rhs2,
int* ridx2,
int rn2,
5969 Rational* rhs3,
int* ridx3,
int rn3 )
6000 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
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]...
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase's dimension to newdim.
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.
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)
The loaded matrix is singular.
DVectorRational diag
Array of pivot elements.
void solveUpdateLeft(Rational *vec)
Temp temp
Temporary storage.
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)
DVectorRational val
hold nonzero values: this is only initialized in the end of the factorization with DEFAULT updates...
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)
DVectorRational s_max
maximum absolute value per row (or -1)
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)
DVectorRational val
values of L vectors
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)
DVectorRational rval
values of rows of L
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
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
DVectorRational val
hold nonzero values
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.
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
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.