32 #define SOPLEX_FACTOR_MARKER 1e-100 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] )
260 for (
int i = 0; i <
thedim; ++i )
283 <<
"LU pivot element is almost zero (< " 285 <<
") - Basis is numerically singular" 314 for ( ring = list->
next; ring != list; ring = ring->
next )
318 if ( l_rbeg[l_row] != n )
324 assert( l_rlen[l_row] <= l_rmax[l_row] );
326 l_rmax[l_row] = l_rlen[l_row];
327 j = i + l_rlen[l_row];
329 for ( ; i < j; ++i, ++n )
332 l_ridx[n] = l_ridx[i];
333 l_rval[n] = l_rval[i];
338 while ( ring != list );
340 goto terminatePackRows;
345 l_rmax[l_row] = l_rlen[l_row];
372 for ( ring = list->
next; ring != list; ring = ring->
next )
376 if ( cbeg[colno] != n )
383 cmax[colno] = clen[colno];
394 while ( ring != list );
396 goto terminatePackColumns;
401 cmax[colno] = clen[colno];
404 terminatePackColumns :
415 assert(
u.
row.
max[p_row] < len );
419 int delta = len -
u.
row.
max[p_row];
424 delta = len -
u.
row.
max[p_row];
434 &&
"ERROR: could not allocate memory for row file" );
458 &&
"ERROR: could not allocate memory for row file" );
475 for ( ; i < k; ++i, ++j )
504 for ( ring = list->
next; ring != list; ring = ring->
next )
508 if ( l_cbeg[l_col] != n )
515 l_cmax[l_col] = l_clen[l_col];
516 j = i + l_clen[l_col];
519 l_cidx[n++] = l_cidx[i];
523 while ( ring != list );
525 goto terminatePackColumns;
530 l_cmax[l_col] = l_clen[l_col];
533 terminatePackColumns :
544 assert(
u.
col.
max[p_col] < len );
548 int delta = len -
u.
col.
max[p_col];
553 delta = len -
u.
col.
max[p_col];
563 &&
"ERROR: could not allocate memory for column file" );
587 &&
"ERROR: could not allocate memory for column file" );
613 assert(
u.
col.
max[p_col] < len );
617 int delta = len -
u.
col.
max[p_col];
622 delta = len -
u.
col.
max[p_col];
630 &&
"ERROR: could not allocate memory for column file" );
653 &&
"ERROR: could not allocate memory for column file" );
712 int i, j, k, h, m, n;
746 for ( i += j - 1; i >= j; --i )
750 h = --( rlen[m] ) + k;
752 while ( ridx[k] != p_col )
777 if ( num > cmax[p_col] )
788 for ( j = 0; j < num; ++j )
796 if (
spxAbs( x ) > l_maxabs )
800 assert( k - cbeg[p_col] < cmax[p_col] );
807 if ( rmax[i] <= rlen[i] )
814 h = rbeg[i] + ( rlen[i] )++;
826 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
843 for ( i = 0; i < dim; ++i )
850 if (
spxAbs( x ) > l_maxabs )
856 clen[p_col] = k - cbeg[p_col];
865 assert( k - cbeg[p_col] < cmax[p_col] );
872 if ( rmax[i] <= rlen[i] )
879 h = rbeg[i] + ( rlen[i] )++;
891 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
893 if ( cbeg[p_col] + cmax[p_col] ==
u.
col.
used )
896 cmax[p_col] = clen[p_col];
912 memmove(&rorig[c], &rorig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
916 for ( i = c; i <= r; ++i )
924 memmove(&corig[c], &corig[c + 1], (
unsigned int)(r - c) *
sizeof(
int));
928 for ( i = c; i <= r; ++i )
940 if ( i < verySparseFactor * ( dim - c ) )
957 for ( i += j - 1; i >= j; --i )
962 m = --( clen[k] ) + cbeg[k];
964 for ( h = m; cidx[h] != rowno; --h )
981 assert(( num == 0 ) || ( nonz != 0 ) );
989 for ( i = 0; i < num; ++i )
990 assert( p_work[corig[nonz[i]]] != 0.0 );
1000 assert( p_work[k] != 0.0 );
1004 x = p_work[k] *
diag[n];
1014 if (
spxAbs( x ) > l_maxabs )
1021 for ( ; j < m; ++j )
1024 Real y = p_work[jj];
1052 diag[rowno] = 1 / x;
1059 if ( rmax[rowno] < num )
1073 for ( i = 0; i < num; ++i )
1086 if (
spxAbs( x ) > l_maxabs )
1097 if ( clen[j] >= cmax[j] )
1104 cval[cbeg[j] + clen[j]] = x;
1106 cidx[cbeg[j] + clen[j]++] = rowno;
1110 rlen[rowno] = n - rbeg[rowno];
1116 for ( i += j - 1; i >= j; --i )
1119 p_work[k] = rval[i];
1120 m = --( clen[k] ) + cbeg[k];
1122 for ( h = m; cidx[h] != rowno; --h )
1139 for ( i = c; i < r; ++i )
1143 if ( p_work[k] != 0.0 )
1146 x = p_work[k] *
diag[n];
1152 if (
spxAbs( x ) > l_maxabs )
1159 for ( ; j < m; ++j )
1160 p_work[ridx[j]] -= x * rval[j];
1183 diag[rowno] = 1 / x;
1193 for ( i = r + 1; i < dim; ++i )
1194 if ( p_work[corig[i]] != 0.0 )
1197 if ( rmax[rowno] < n )
1211 for ( i = r + 1; i < dim; ++i )
1218 if (
spxAbs( x ) > l_maxabs )
1229 if ( clen[j] >= cmax[j] )
1236 cval[cbeg[j] + clen[j]] = x;
1238 cidx[cbeg[j] + clen[j]++] = rowno;
1242 rlen[rowno] = n - rbeg[rowno];
1253 i = rbeg[rowno] + --( rlen[rowno] );
1254 diag[rowno] = 1 / rval[i];
1256 for ( j = i = --( clen[p_col] ) + cbeg[p_col]; cidx[i] != rowno; --i )
1283 assert( p_work[p_col] != 0.0 );
1284 rezi = 1 / p_work[p_col];
1285 p_work[p_col] = 0.0;
1292 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1295 lval[ll] = rezi * p_work[j];
1302 lval[ll] = 1 - rezi;
1305 for ( --i; i >= 0; --i )
1309 lval[ll] = x = rezi * p_work[j];
1331 assert( p_work[p_col] != 0.0 );
1332 rezi = 1 / p_work[p_col];
1338 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1341 lval[ll] = rezi * p_work[j];
1347 lval[ll] = 1 - rezi;
1350 for ( --i; i >= 0; --i )
1354 lval[ll] = x = rezi * p_work[j];
1396 Dring *rring, *lastrring;
1397 Dring *cring, *lastcring;
1407 for (
int i = 0; i <
thedim; i++ )
1412 for (
int i = 0; i <
thedim; i++ )
1423 for (
int j = 0; j < k; ++j )
1456 lastrring->
next = rring;
1464 lastcring->
next = cring;
1468 for (
int i = 0; i <
thedim; i++ )
1474 rring->
prev = lastrring;
1475 lastrring->
next = rring;
1480 cring->
prev = lastcring;
1481 lastcring->
next = cring;
1505 for (
int i = 0; i <
thedim; i++ )
1515 for (
int j = 0; j < psv->
size() && nnonzeros <= 1; j++ )
1522 if ( nnonzeros == 0 )
1530 if ( nnonzeros == 1 )
1536 for ( j = 0;
isZero( psv->
value( j ), eps ); j++ )
1539 assert( j < psv->size() );
1549 x = psv->
value( j );
1572 assert( nnonzeros >= 2 );
1575 for (
int j = 0; j < psv->
size(); j++ )
1577 x = psv->
value( j );
1582 k = psv->
index( j );
1601 assert( nnonzeros >= 2 );
1622 int p_col, p_row, newrow;
1638 assert( p_row >= 0 );
1642 for ( j = 0; j < len; ++j )
1651 for ( k = n;
u.
col.
idx[k] != p_row; ++k )
1669 if ( rperm[newrow] >= 0 )
1679 for ( k = n;
u.
row.
idx[k] != p_col; --k )
1715 int p_row, p_col, len, rs, lk;
1724 for ( i = 0; i <
thedim; ++i )
1726 if ( rperm[i] < 0 &&
u.
row.
len[i] == 1 )
1742 setPivot( rs, p_col, p_row, pval );
1752 i = (
u.
col.
len[p_col] -= i );
1754 for ( ; i < len; ++i )
1765 for ( j = k;
u.
row.
idx[j] != p_col; --j )
1825 for ( i = 0; i <
thedim; ++i )
1910 for ( ; ( r = idx[i] ) != prow; ++i )
1917 for ( j = k;
u.
row.
idx[j] != pcol; --j )
1948 assert( i < len &&
"ERROR: pivot column does not contain pivot row" );
1950 for ( ++i; i < len; ++i )
1958 for ( j = k;
u.
row.
idx[j] != pcol; --j )
2020 for ( i = j; ( c =
u.
row.
idx[i] ) != pcol; --i )
2024 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2050 for ( --i; i >= j; --i )
2055 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2100 if ( candidates > 4 )
2115 len =
u.
row.
len[rw] + beg - 1;
2125 for ( i = len - 1; i >= beg; --i )
2132 l_maxabs *= threshold;
2138 for ( i = len; i >= beg; --i )
2144 if ( j < mkwtz &&
spxAbs( x ) > l_maxabs )
2160 len =
u.
col.
len[cl] + beg - 1;
2168 for ( i = len; i >= beg; --i )
2201 for ( kk = m +
u.
row.
len[k] - 1; kk >= m; --kk )
2214 for ( --kk; kk > m; --kk )
2223 l_maxabs *= threshold;
2225 if (
spxAbs( x ) > l_maxabs )
2231 if ( j <= count + 1 )
2260 if ( pr->
idx == rw || pr->
mkwtz >= mkwtz )
2266 if ( pr->
idx != rw )
2274 if ( num >= candidates )
2307 int c, i, j, k, ll, m, n;
2316 for ( j = m;
u.
row.
idx[j] != pcol; --j )
2337 for ( j = m - 1; j >= n; --j )
2366 for ( i = k;
u.
col.
idx[i] != r; --i )
2379 if ( ll + fill >
u.
row.
max[r] )
2436 int i, j, k, m = -1;
2441 int plen = --(
u.
row.
len[prow] );
2442 int pend = pbeg + plen;
2471 for ( i = pbeg; i < pend; ++i )
2479 for ( k = m;
u.
col.
idx[k] != prow; ++k )
2497 updateRow( m, lv++, prow, pcol, pval, eps );
2504 for ( ++i; i < m; ++i )
2516 for ( i =
u.
row.
start[prow], pend = i + plen; i < pend; ++i )
2531 const Real threshold )
2551 for ( i = 0; i <
thedim; ++i )
2570 "ERROR: no pivot element selected" );
2573 pivot = pivot->
next )
2593 "ERROR: one row must be left" );
2595 "ERROR: one col must be left" );
2617 for ( i = 0; i <
thedim; i++ )
2622 for ( i = 0; i <
thedim; i++ )
2633 assert(( *idx >= 0 ) && ( *idx < thedim ) );
2637 assert(( k >= 0 ) && ( k <
u.
col.
size ) );
2723 for ( i =
thedim; i--; *l_rbeg++ = 0 )
2725 *rorig++ = *rrorig++;
2726 *rperm++ = *rrperm++;
2731 l_rbeg =
l.
rbeg + 1;
2733 for ( i = mem; i--; )
2738 for ( m = 0, i =
thedim; i--; l_rbeg++ )
2747 l_rbeg =
l.
rbeg + 1;
2749 for ( i = j = 0; i < vecs; ++i )
2754 for ( ; j < beg[i + 1]; j++ )
2756 k = l_rbeg[*idx++]++;
2765 assert(
l.
rbeg[0] == 0 );
2838 for ( i = 0; i <
thedim; ++i )
2841 std::cout <<
"DCLUFA01 diag[" << i <<
"]: [" <<
col.
orig[
row.
perm[i]]
2842 <<
"] = " <<
diag[i] << std::endl;
2844 for ( j = 0; j <
u.
row.
len[i]; ++j )
2845 std::cout <<
"DCLUFA02 u[" << i <<
"]: [" 2852 for ( i = 0; i <
thedim; ++i )
2857 std::cout <<
"DCLUFA03 l[" << i <<
"]" << std::endl;
2860 std::cout <<
"DCLUFA04 l[" << k -
l.
start[j]
2861 <<
"]: [" <<
l.
idx[k]
2862 <<
"] = " <<
l.
val[k] << std::endl;
2914 if ( size >
l.
size )
2932 int* p_lrow =
l.
row;
2937 assert( p_len > 0 &&
"ERROR: no empty columns allowed in L vectors" );
2953 #ifdef ENABLE_CONSISTENCY_CHECKS 2970 assert( ring->
idx >= 0 );
2972 assert( ring->
prev->
next == ring );
2995 assert( ring->
idx >= 0 );
2997 assert( ring->
prev->
next == ring );
3012 for ( i = 0; i <
thedim; ++i )
3021 for ( i = 0; i <
thedim; ++i )
3048 for ( i = 0; i <
thedim; ++i )
3070 for ( i = 0; i < thedim -
temp.
stage; ++i )
3084 #if 0 // Stimmt leider nicht 3089 for ( i = 0; i <
thedim; i++ )
3091 assert(
l.
rbeg[i] >= 0 );
3095 assert(
l.
rorig[i] >= 0 );
3096 assert(
l.
rorig[i] < thedim );
3097 assert(
l.
rperm[i] >= 0 );
3098 assert(
l.
rperm[i] < thedim );
3100 assert(
l.
ridx[i] >= 0 );
3101 assert(
l.
ridx[i] < thedim );
3106 #endif // CONSISTENCY_CHECKS 3114 for (
int i =
thedim - 1; i >= 0; i-- )
3118 Real x = wrk[c] =
diag[r] * vec[r];
3135 int *cidx, *clen, *cbeg;
3152 for ( i =
thedim - 1; i >= 0; --i )
3155 x =
diag[r] * rhs[r];
3162 val = &cval[cbeg[c]];
3163 idx = &cidx[cbeg[c]];
3167 rhs[*idx++] -= x * ( *val++ );
3179 int *cidx, *clen, *cbeg;
3194 for ( i =
thedim - 1; i >= 0; --i )
3198 p_work1[c] = x1 =
diag[r] * vec1[r];
3199 p_work2[c] = x2 =
diag[r] * vec2[r];
3200 vec1[r] = vec2[r] = 0;
3202 if ( x1 != 0.0 && x2 != 0.0 )
3204 val = &cval[cbeg[c]];
3205 idx = &cidx[cbeg[c]];
3210 vec1[*idx] -= x1 * ( *val );
3211 vec2[*idx++] -= x2 * ( *val++ );
3217 val = &cval[cbeg[c]];
3218 idx = &cidx[cbeg[c]];
3222 vec1[*idx++] -= x1 * ( *val++ );
3227 val = &cval[cbeg[c]];
3228 idx = &cidx[cbeg[c]];
3232 vec2[*idx++] -= x2 * ( *val++ );
3239 int* nonz,
Real eps )
3243 int *cidx, *clen, *cbeg;
3244 bool notzero1, notzero2;
3261 for ( i =
thedim - 1; i >= 0; --i )
3265 p_work1[c] = x1 =
diag[r] * vec1[r];
3266 p_work2[c] = x2 =
diag[r] * vec2[r];
3267 vec1[r] = vec2[r] = 0;
3271 if ( notzero1 && notzero2 )
3275 val = &cval[cbeg[c]];
3276 idx = &cidx[cbeg[c]];
3281 vec1[*idx] -= x1 * ( *val );
3282 vec2[*idx++] -= x2 * ( *val++ );
3291 val = &cval[cbeg[c]];
3292 idx = &cidx[cbeg[c]];
3296 vec1[*idx++] -= x1 * ( *val++ );
3302 val = &cval[cbeg[c]];
3303 idx = &cidx[cbeg[c]];
3307 vec2[*idx++] -= x2 * ( *val++ );
3325 int *lrow, *lidx, *idx;
3335 for ( i = 0; i < end; ++i )
3337 if (( x = vec[lrow[i]] ) != 0.0 )
3339 MSG_DEBUG( std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl; )
3345 for ( j = lbeg[i + 1]; j > k; --j )
3347 MSG_DEBUG( std::cout <<
" -> y" << *idx <<
" -= " << x <<
" * " << *val <<
" = " << x * ( *val ) <<
" -> " << vec[*idx] - x * ( *val ) << std::endl; )
3348 vec[*idx++] -= x * ( *val++ );
3355 MSG_DEBUG( std::cout <<
"performing FT updates..." << std::endl; )
3359 for ( ; i < end; ++i )
3366 for ( j = lbeg[i + 1]; j > k; --j )
3367 x += vec[*idx++] * ( *val++ );
3371 MSG_DEBUG( std::cout <<
"y" << lrow[i] <<
"=" << vec[lrow[i]] << std::endl; )
3374 MSG_DEBUG( std::cout <<
"finished FT updates." << std::endl; )
3385 int *lrow, *lidx, *idx;
3395 for ( i = 0; i < end; ++i )
3400 if ( x1 != 0 && x2 != 0.0 )
3406 for ( j = lbeg[i + 1]; j > k; --j )
3408 vec1[*idx] -= x1 * ( *val );
3409 vec2[*idx++] -= x2 * ( *val++ );
3419 for ( j = lbeg[i + 1]; j > k; --j )
3420 vec1[*idx++] -= x1 * ( *val++ );
3429 for ( j = lbeg[i + 1]; j > k; --j )
3430 vec2[*idx++] -= x2 * ( *val++ );
3438 for ( ; i < end; ++i )
3446 for ( j = lbeg[i + 1]; j > k; --j )
3448 x1 += vec1[*idx] * ( *val );
3449 x2 += vec2[*idx++] * ( *val++ );
3452 vec1[lrow[i]] -= x1;
3454 vec2[lrow[i]] -= x2;
3465 int *lrow, *lidx, *idx;
3479 if (( x = vec[lrow[i]] ) != 0.0 )
3485 for ( j = lbeg[i + 1]; j > k; --j )
3486 vec[*idx++] -= x * ( *val++ );
3517 if ( x1 != 0.0 && x2 != 0.0 )
3523 for ( j = lbeg[i + 1]; j > k; --j )
3525 vec1[*idx] -= x1 * ( *val );
3526 vec2[*idx++] -= x2 * ( *val++ );
3536 for ( j = lbeg[i + 1]; j > k; --j )
3537 vec1[*idx++] -= x1 * ( *val++ );
3546 for ( j = lbeg[i + 1]; j > k; --j )
3547 vec2[*idx++] -= x2 * ( *val++ );
3553 Real* rhs,
Real* forest,
int* forestNum,
int* forestIdx )
3561 for (
int i = 0; i <
thedim; i++ )
3565 n += rhs[i] != 0.0 ? 1 : 0;
3606 for (
int i = 0; i <
thedim; i++ )
3609 forest[i] = rhs1[i];
3610 n += rhs1[i] != 0.0 ? 1 : 0;
3650 int *ridx, *rlen, *rbeg, *idx;
3661 for ( i = 0; i <
thedim; ++i )
3676 for (
int m = rlen[r]; m != 0; --m )
3677 vec[*idx++] -= x * ( *val++ );
3687 for (
int i = 0; i <
thedim; ++i )
3692 for (
int m =
u.
row.
start[r]; m < end; m++ )
3694 std::cout <<
u.
row.
val[m] <<
" ";
3701 for (
int i = 0; i <
thedim; ++i )
3714 #if defined(WITH_WARNINGS) || !defined(NDEBUG) 3731 for (
int m =
u.
row.
start[r]; m < end; m++ )
3753 int *ridx, *rlen, *rbeg, *idx;
3764 for ( i = 0; i <
thedim; ++i )
3771 if (( x1 != 0.0 ) && ( x2 != 0.0 ) )
3781 for (
int m = rlen[r]; m != 0; --m )
3783 vec1[*idx] -= x1 * ( *val );
3784 vec2[*idx] -= x2 * ( *val++ );
3797 for (
int m = rlen[r]; m != 0; --m )
3798 vec1[*idx++] -= x1 * ( *val++ );
3809 for (
int m = rlen[r]; m != 0; --m )
3810 vec2[*idx++] -= x2 * ( *val++ );
3827 int *lidx, *idx, *lrow;
3851 for ( j = lbeg[i + 1]; j > k; --j )
3853 vec1[*idx] -= x1 * ( *val );
3854 vec2[*idx++] -= x2 * ( *val++ );
3863 for ( j = lbeg[i + 1]; j > k; --j )
3864 vec1[*idx++] -= x1 * ( *val++ );
3874 for ( j = lbeg[i + 1]; j > k; --j )
3875 vec2[*idx++] -= x2 * ( *val++ );
3910 for ( ; i >= 0; --i )
3918 for ( j = lbeg[i + 1]; j > k; --j )
3920 x1 += vec1[*idx] * ( *val );
3921 x2 += vec2[*idx++] * ( *val++ );
3924 vec1[lrow[i]] -= x1;
3926 vec2[lrow[i]] -= x2;
3935 x1not0 = ( x1 != 0 );
3936 x2not0 = ( x2 != 0 );
3938 if ( x1not0 && x2not0 )
3941 j = rbeg[r + 1] - k;
3948 vec1[*idx] -= x1 * *val;
3949 vec2[*idx++] -= x2 * *val++;
3956 j = rbeg[r + 1] - k;
3963 vec1[*idx++] -= x1 * *val++;
3970 j = rbeg[r + 1] - k;
3977 vec2[*idx++] -= x2 * *val++;
3990 int *idx, *lidx, *lrow, *lbeg;
4001 if (( x = vec[lrow[i]] ) != 0.0 )
4007 for ( j = lbeg[i + 1]; j > k; --j )
4008 vec[*idx++] -= x * ( *val++ );
4033 for (
int j = lbeg[i + 1]; j > k; --j )
4034 x += vec[*idx++] * ( *val++ );
4040 for (
int i =
thedim - 1; i >= 0; --i )
4047 for (
int k =
l.
rbeg[r]; k <
l.
rbeg[r + 1]; k++ )
4051 assert(
l.
rperm[j] < i );
4053 vec[j] -= x *
l.
rval[k];
4058 #endif // WITH_L_ROWS 4089 for ( j = lbeg[i + 1]; j > k; --j )
4090 x += vec[*idx++] * ( *val++ );
4106 j = rbeg[r + 1] - k;
4113 vec[*idx++] -= x * *val++;
4130 int *lrow, *lidx, *idx;
4149 for ( j = lbeg[i + 1]; j > k; --j )
4150 x += vec[*idx++] * ( *val++ );
4161 int *lrow, *lidx, *idx;
4181 for ( j = lbeg[i + 1]; j > k; --j )
4183 x1 += vec1[*idx] * ( *val );
4184 x2 += vec2[*idx++] * ( *val++ );
4187 vec1[lrow[i]] -= x1;
4189 vec2[lrow[i]] -= x2;
4198 int *lrow, *lidx, *idx;
4213 assert( k >= 0 && k <
l.
size );
4218 for ( j = lbeg[i + 1]; j > k; --j )
4220 assert( *idx >= 0 && *idx <
thedim );
4221 x += vec[*idx++] * ( *val++ );
4309 Real* vec,
int* vecidx,
4310 Real* rhs,
int* rhsidx,
int rhsn )
4313 int i, j, k, n, r, c;
4314 int *rorig, *corig, *cperm;
4315 int *ridx, *rlen, *rbeg, *idx;
4325 for ( i = 0; i < rhsn; )
4341 assert( i >= 0 && i <
thedim );
4343 assert( c >= 0 && c <
thedim );
4350 assert( r >= 0 && r <
thedim );
4355 assert( k >= 0 && k <
u.
row.
size );
4359 for (
int m = rlen[r]; m; --m )
4362 assert( j >= 0 && j <
thedim );
4367 y = -x * ( *val++ );
4377 y -= x * ( *val++ );
4389 Real* rhs,
int* rhsidx,
int rhsn )
4393 int *rorig, *corig, *cperm;
4394 int *ridx, *rlen, *rbeg, *idx;
4404 for ( i = 0; i < rhsn; )
4418 assert( i >= 0 && i <
thedim );
4420 assert( c >= 0 && c <
thedim );
4427 assert( r >= 0 && r <
thedim );
4431 assert( k >= 0 && k <
u.
row.
size );
4435 for (
int m = rlen[r]; m; --m )
4438 assert( j >= 0 && j <
thedim );
4443 y = -x * ( *val++ );
4453 y -= x * ( *val++ );
4467 int *idx, *lidx, *lrow, *lbeg;
4477 assert( i >= 0 && i <
l.
size );
4479 if (( x = vec[lrow[i]] ) != 0.0 )
4482 assert( k >= 0 && k <
l.
size );
4486 for ( j = lbeg[i + 1]; j > k; --j )
4489 assert( m >= 0 && m <
thedim );
4494 y = -x * ( *val++ );
4504 y -= x * ( *val++ );
4520 int *idx, *lidx, *lrow, *lbeg;
4530 if (( x = vec[lrow[i]] ) != 0.0 )
4532 assert( i >= 0 && i <
l.
size );
4534 assert( k >= 0 && k <
l.
size );
4538 for ( j = lbeg[i + 1]; j > k; --j )
4540 assert( *idx >= 0 && *idx <
thedim );
4541 vec[*idx++] -= x * ( *val++ );
4568 #pragma warn "Not yet implemented, define WITH_L_ROWS" 4574 for ( ; i >= 0; --i )
4581 for ( j = lbeg[i + 1]; j > k; --j )
4582 x += vec[*idx++] * ( *val++ );
4590 for ( i = 0; i < rn; )
4606 j = rbeg[r + 1] - k;
4612 assert(
l.
rperm[*idx] < i );
4637 for ( i = 0; i < n; ++i )
4670 for ( ; i >= 0; --i )
4673 assert( k >= 0 && k <
l.
size );
4678 for ( j = lbeg[i + 1]; j > k; --j )
4680 assert( *idx >= 0 && *idx <
thedim );
4681 x += vec[*idx++] * ( *val++ );
4696 j = rbeg[r + 1] - k;
4702 assert(
l.
rperm[*idx] < i );
4703 vec[*idx++] -= x * *val++;
4734 int *lrow, *lidx, *idx;
4745 for( i = 0; i < end; ++i )
4757 for( j = lbeg[i + 1]; j > k; --j )
4759 assert(*idx >= 0 && *idx <
thedim);
4771 for( ; i < end; ++i )
4778 for( j = lbeg[i + 1]; j > k; --j )
4780 assert(*idx >= 0 && *idx <
thedim);
4781 x += vec[*idx++] * ( *val++ );
4795 Real* vec,
int* ridx,
int& rn,
Real eps,
4796 Real* vec2,
int* ridx2,
int& rn2,
Real eps2 )
4802 int *lrow, *lidx, *idx;
4813 for ( i = 0; i < end; ++i )
4829 for( j = lbeg[i + 1]; j > k; --j )
4831 assert(*idx >= 0 && *idx <
thedim);
4841 for( j = lbeg[i + 1]; j > k; --j )
4843 assert(*idx >= 0 && *idx <
thedim);
4857 for( j = lbeg[i + 1]; j > k; --j )
4859 assert(*idx >= 0 && *idx <
thedim);
4871 for( ; i < end; ++i )
4878 for( j = lbeg[i + 1]; j > k; --j )
4880 assert( *idx >= 0 && *idx <
thedim );
4881 x += vec[*idx] * ( *val );
4882 x2 += vec2[*idx++] * ( *val++ );
4899 Real* vec,
int* ridx,
int& rn,
Real eps,
4900 Real* vec2,
int* ridx2,
int& rn2,
Real eps2,
4901 Real* vec3,
int* ridx3,
int& rn3,
Real eps3 )
4907 int *lrow, *lidx, *idx;
4917 for( i = 0; i < end; ++i )
4935 for( j = lbeg[i + 1]; j > k; --j )
4937 assert( *idx >= 0 && *idx <
thedim );
4948 for( j = lbeg[i + 1]; j > k; --j )
4950 assert( *idx >= 0 && *idx <
thedim );
4961 for ( j = lbeg[i + 1]; j > k; --j )
4963 assert( *idx >= 0 && *idx <
thedim );
4973 for ( j = lbeg[i + 1]; j > k; --j )
4975 assert( *idx >= 0 && *idx <
thedim );
4991 for( j = lbeg[i + 1]; j > k; --j )
4993 assert( *idx >= 0 && *idx <
thedim );
5003 for( j = lbeg[i + 1]; j > k; --j )
5005 assert( *idx >= 0 && *idx <
thedim );
5019 for ( j = lbeg[i + 1]; j > k; --j )
5021 assert( *idx >= 0 && *idx <
thedim );
5033 for ( ; i < end; ++i )
5040 for ( j = lbeg[i + 1]; j > k; --j )
5042 assert( *idx >= 0 && *idx <
thedim );
5043 x += vec[*idx] * ( *val );
5044 x2 += vec2[*idx] * ( *val );
5045 x3 += vec3[*idx++] * ( *val++ );
5063 Real* rhs,
int* ridx,
int rn,
Real eps )
5065 int i, j, k, r, c, n;
5068 int *cidx, *clen, *cbeg;
5091 assert( i >= 0 && i <
thedim );
5093 assert( r >= 0 && r <
thedim );
5095 x =
diag[r] * rhs[r];
5101 assert( c >= 0 && c <
thedim );
5104 val = &cval[cbeg[c]];
5105 idx = &cidx[cbeg[c]];
5110 assert( *idx >= 0 && *idx <
thedim );
5112 assert( k >= 0 && k <
thedim );
5117 y = -x * ( *val++ );
5127 y -= x * ( *val++ );
5133 if ( rn > i*verySparseFactor4right )
5135 for ( i = *ridx; i >= 0; --i )
5138 assert( r >= 0 && r <
thedim );
5139 x =
diag[r] * rhs[r];
5145 assert( c >= 0 && c <
thedim );
5148 val = &cval[cbeg[c]];
5149 idx = &cidx[cbeg[c]];
5154 assert( *idx >= 0 && *idx <
thedim );
5155 rhs[*idx++] -= x * ( *val++ );
5170 Real* rhs,
int* ridx,
int rn,
Real eps )
5175 int *cidx, *clen, *cbeg;
5193 if ( rn > *ridx * verySparseFactor4right )
5195 for ( i = *ridx; i >= 0; --i )
5197 assert( i >= 0 && i <
thedim );
5199 assert( r >= 0 && r <
thedim );
5200 x =
diag[r] * rhs[r];
5207 val = &cval[cbeg[c]];
5208 idx = &cidx[cbeg[c]];
5213 assert( *idx >= 0 && *idx <
thedim );
5214 rhs[*idx++] -= x * ( *val++ );
5226 assert( i >= 0 && i <
thedim );
5230 assert( r >= 0 && r <
thedim );
5232 x =
diag[r] * rhs[r];
5240 val = &cval[cbeg[c]];
5241 idx = &cidx[cbeg[c]];
5247 assert( k >= 0 && k <
thedim );
5252 y = -x * ( *val++ );
5262 y -= x * ( *val++ );
5273 Real* vec,
int* vidx,
Real* rhs,
int* ridx,
int rn,
Real eps,
5274 Real* vec2,
Real* rhs2,
int* ridx2,
int rn2,
Real eps2 )
5276 int i, j, k, r, c, n;
5279 int *cidx, *clen, *cbeg;
5298 while ( rn + rn2 > 0 )
5308 if ( *ridx2 > *ridx )
5311 if ( *ridx2 < *ridx )
5319 assert( i >= 0 && i <
thedim );
5322 assert( r >= 0 && r <
thedim );
5324 x =
diag[r] * rhs[r];
5325 x2 =
diag[r] * rhs2[r];
5335 val = &cval[cbeg[c]];
5336 idx = &cidx[cbeg[c]];
5344 assert( k >= 0 && k <
thedim );
5349 y2 = -x2 * ( *val );
5359 y2 -= x2 * ( *val );
5367 y = -x * ( *val++ );
5377 y -= x * ( *val++ );
5388 assert( k >= 0 && k <
thedim );
5393 y = -x * ( *val++ );
5403 y -= x * ( *val++ );
5414 assert( c >= 0 && c <
thedim );
5416 val = &cval[cbeg[c]];
5417 idx = &cidx[cbeg[c]];
5423 assert( k >= 0 && k <
thedim );
5428 y2 = -x2 * ( *val++ );
5438 y2 -= x2 * ( *val++ );
5444 if ( rn + rn2 > i*verySparseFactor4right )
5446 if ( *ridx > *ridx2 )
5451 for ( ; i >= 0; --i )
5455 assert( r >= 0 && r <
thedim );
5456 x =
diag[r] * rhs[r];
5457 x2 =
diag[r] * rhs2[r];
5464 assert( c >= 0 && c <
thedim );
5466 val = &cval[cbeg[c]];
5467 idx = &cidx[cbeg[c]];
5477 assert( *idx >= 0 && *idx <
thedim );
5478 rhs[*idx] -= x * ( *val );
5479 rhs2[*idx++] -= x2 * ( *val++ );
5486 assert( *idx >= 0 && *idx <
thedim );
5487 rhs2[*idx++] -= x2 * ( *val++ );
5495 assert( c >= 0 && c <
thedim );
5498 val = &cval[cbeg[c]];
5499 idx = &cidx[cbeg[c]];
5504 assert( *idx >= 0 && *idx <
thedim );
5505 rhs[*idx++] -= x * ( *val++ );
5523 int *lrow, *lidx, *idx;
5536 assert( i >= 0 && i <
thedim );
5542 assert( k >= 0 && k <
l.
size );
5546 for ( j = lbeg[i + 1]; j > k; --j )
5548 int m = ridx[n] = *idx++;
5549 assert( m >= 0 && m <
thedim );
5551 n += ( y == 0 ) ? 1 : 0;
5552 y = y - x * ( *val++ );
5567 int *lrow, *lidx, *idx;
5580 assert( i >= 0 && i <
thedim );
5582 if (( x = vec[lrow[i]] ) != 0.0 )
5585 assert( k >= 0 && k <
l.
size );
5589 for ( j = lbeg[i + 1]; j > k; --j )
5591 assert( *idx >= 0 && *idx <
thedim );
5592 vec[*idx++] -= x * ( *val++ );
5600 Real* vec,
int* idx,
5601 Real* rhs,
int* ridx,
int rn,
5602 Real* forest,
int* forestNum,
int* forestIdx )
5605 assert(rn >= 0 && rn <=
thedim);
5615 int* it = forestIdx;
5619 for ( i = j = 0; i < rn; ++i )
5622 assert( k >= 0 && k <
thedim );
5634 *forestNum = rn = j;
5644 for ( i = j = 0; i < rn; ++i )
5647 assert( k >= 0 && k <
thedim );
5668 Real* vec,
int* idx,
5669 Real* rhs,
int* ridx,
int rn,
5671 Real* rhs2,
int* ridx2,
int rn2,
5672 Real* forest,
int* forestNum,
int* forestIdx )
5675 assert(rn >= 0 && rn <=
thedim);
5676 assert(rn2 >= 0 && rn2 <=
thedim);
5686 int* it = forestIdx;
5690 for ( i = j = 0; i < rn; ++i )
5693 assert( k >= 0 && k <
thedim );
5705 *forestNum = rn = j;
5715 for ( i = j = 0; i < rn; ++i )
5718 assert( k >= 0 && k <
thedim );
5730 if ( rn2 >
thedim*verySparseFactor4right )
5745 for ( i = j = 0; i < rn2; ++i )
5748 assert( k >= 0 && k <
thedim );
5789 Real* rhs,
int* ridx,
int& rn,
5791 Real* rhs2,
int* ridx2,
int& rn2,
5792 Real* forest,
int* forestNum,
int* forestIdx )
5796 assert(rn >= 0 && rn <=
thedim);
5797 assert(rn2 >= 0 && rn2 <=
thedim);
5806 int* it = forestIdx;
5808 for ( i = j = 0; i < rn; ++i )
5811 assert( k >= 0 && k <
thedim );
5823 *forestNum = rn = j;
5827 for ( i = j = 0; i < rn; ++i )
5830 assert( k >= 0 && k <
thedim );
5842 for ( i = j = 0; i < rn2; ++i )
5845 assert( k >= 0 && k <
thedim );
5858 rn2 =
vSolveUright( vec2, idx2, rhs2, ridx2, rn2, eps2 );
5869 Real* vec,
int* idx,
5870 Real* rhs,
int* ridx,
int rn,
5872 Real* rhs2,
int* ridx2,
int rn2,
5874 Real* rhs3,
int* ridx3,
int rn3,
5875 Real* forest,
int* forestNum,
int* forestIdx )
5878 vSolveLright3(rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2, rhs3, ridx3, rn3, eps3);
5879 assert(rn >= 0 && rn <=
thedim);
5880 assert(rn2 >= 0 && rn2 <=
thedim);
5881 assert(rn3 >= 0 && rn3 <=
thedim);
5891 int* it = forestIdx;
5895 for ( i = j = 0; i < rn; ++i )
5898 assert( k >= 0 && k <
thedim );
5910 *forestNum = rn = j;
5920 for ( i = j = 0; i < rn; ++i )
5923 assert( k >= 0 && k <
thedim );
5935 if ( rn2 >
thedim*verySparseFactor4right )
5947 for ( i = j = 0; i < rn2; ++i )
5950 assert( k >= 0 && k <
thedim );
5969 if ( rn3 >
thedim*verySparseFactor4right )
5981 for ( i = j = 0; i < rn3; ++i )
5984 assert( k >= 0 && k <
thedim );
6019 Real* rhs,
int* ridx,
int& rn,
6021 Real* rhs2,
int* ridx2,
int& rn2,
6023 Real* rhs3,
int* ridx3,
int& rn3,
6024 Real* forest,
int* forestNum,
int* forestIdx )
6026 vSolveLright3( rhs, ridx, rn, eps, rhs2, ridx2, rn2, eps2, rhs3, ridx3, rn3, eps3 );
6027 assert( rn >= 0 && rn <=
thedim );
6028 assert( rn2 >= 0 && rn2 <=
thedim );
6029 assert( rn3 >= 0 && rn3 <=
thedim );
6038 int* it = forestIdx;
6040 for ( i = j = 0; i < rn; ++i )
6043 assert( k >= 0 && k <
thedim );
6055 *forestNum = rn = j;
6059 for ( i = j = 0; i < rn; ++i )
6062 assert( k >= 0 && k <
thedim );
6074 for ( i = j = 0; i < rn2; ++i )
6077 assert( k >= 0 && k <
thedim );
6088 for ( i = j = 0; i < rn3; ++i )
6091 assert( k >= 0 && k <
thedim );
6103 rn2 =
vSolveUright( vec2, idx2, rhs2, ridx2, rn2, eps2 );
6104 rn3 =
vSolveUright( vec3, idx3, rhs3, ridx3, rn3, eps3 );
6116 Real* rhs,
int* ridx,
int rn )
6119 assert( rn >= 0 && rn <=
thedim );
6121 if ( rn >
thedim*verySparseFactor4right )
6135 for ( i = j = 0; i < rn; ++i )
6138 assert( k >= 0 && k <
thedim );
6168 Real* vec,
int* idx,
6169 Real* rhs,
int* ridx,
int rn )
6175 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6179 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6196 Real* vec,
int* idx,
6197 Real* rhs,
int* ridx,
int rn,
6199 Real* rhs2,
int* ridx2,
int rn2 )
6205 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6211 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6225 Real* vec,
int* idx,
6226 Real* rhs,
int* ridx,
int& rn,
6227 Real* vec2,
int* idx2,
6228 Real* rhs2,
int* ridx2,
int& rn2 )
6233 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6235 rn2 =
solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6239 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6241 rn2 =
solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6252 Real* vec,
int* idx,
6253 Real* rhs,
int* ridx,
int rn,
6255 Real* rhs2,
int* ridx2,
int rn2,
6257 Real* rhs3,
int* ridx3,
int rn3 )
6263 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6271 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6288 Real* vec,
int* idx,
6289 Real* rhs,
int* ridx,
int& rn,
6290 Real* vec2,
int* idx2,
6291 Real* rhs2,
int* ridx2,
int& rn2,
6292 Real* vec3,
int* idx3,
6293 Real* rhs3,
int* ridx3,
int& rn3 )
6298 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6300 rn2 =
solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6302 rn3 =
solveUleft( eps, vec3, idx3, rhs3, ridx3, rn3 );
6306 rn =
solveUleft( eps, vec, idx, rhs, ridx, rn );
6308 rn2 =
solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6310 rn3 =
solveUleft( eps, vec3, idx3, rhs3, ridx3, rn3 );
6322 Real* rhs2,
int* ridx2,
int rn2 )
void solveUleft(Real *work, Real *vec)
Rational spxAbs(const Rational &r)
Absolute.
int * len
used nonzeros per row vectors
bool isNotZero(Real a, Real eps=Param::epsilon())
returns true iff |a| > eps
void solveLleft(Real *vec) const
int updateType
type of updates to be used.
void solveUpdateRight2(Real *vec1, Real *vec2)
void updateNoClear(int p_col, const Real *p_work, const int *p_idx, int num)
int solveUright2eps(Real *work1, Real *vec1, Real *work2, Real *vec2, int *nonz, Real eps)
void vSolveLright3(Real *vec, int *ridx, int &rn, Real eps, Real *vec2, int *ridx2, int &rn2, Real eps2, Real *vec3, int *ridx3, int &rn3, Real eps3)
int solveRight2update(Real *vec1, Real *vec2, Real *rhs1, Real *rhs2, int *nonz, Real eps, Real *forest, int *forestNum, int *forestIdx)
Dring list
Double linked ringlist of vector indices in the order they appear in the column file.
int * max
maximum available nonzeros per row: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
static const Real verySparseFactor4left
int vSolveRight4update(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *forest, int *forestNum, int *forestIdx)
THREADLOCAL const Real infinity
void forestMinColMem(int size)
Memory allocation routines.
void vSolveRight4update3sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real eps2, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real eps3, Real *vec3, int *idx3, Real *rhs3, int *ridx3, int &rn3, Real *forest, int *forestNum, int *forestIdx)
sparse version of above method
Pring * pivot_rowNZ
lists for rows to number of nonzeros
int size() const
Number of used indices.
struct soplex::CLUFactor::U::Col col
int * start
starting positions in val and idx
int solveLleftForest(Real *vec, int *, Real)
Real * rval
values of rows of L
Exception classes for SoPlex.
#define init2DR(elem, ring)
void solveUleftNoNZ(Real eps, Real *vec, Real *rhs, int *rhsidx, int rhsn)
Real maxabs
maximum abs number in L and U
int vSolveRight4update3(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real eps2, Real *rhs2, int *ridx2, int rn2, Real *vec3, Real eps3, Real *rhs3, int *ridx3, int rn3, Real *forest, int *forestNum, int *forestIdx)
void vSolveLeft3sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real *vec3, int *idx3, Real *rhs3, int *ridx3, int &rn3)
sparse version of solving 3 systems of equations
Real * s_max
maximum absolute value per row (or -1)
static void enQueueMin(int *heap, int *size, int elem)
int size
size of array idx
R & value(int n)
Reference to value of n 'th nonzero.
int * rorig
original row permutation
void update(int p_col, Real *p_work, const int *p_idx, int num)
void solveUright(Real *wrk, Real *vec) const
Dring * elem
Array of ring elements.
void solveLeft(Real *vec, Real *rhs)
#define ASSERT_WARN(prefix, expr)
Macro to turn some assertions into warnings.
Real initMaxabs
maximum abs number in initail Matrix
int mkwtz
markowitz number of pivot
Pring * pivot_col
column index handlers for Real linked list
int vSolveLeft2(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real *rhs2, int *ridx2, int rn2)
Real * diag
Array of pivot elements.
Perm row
row permutation matrices
Dring * elem
Array of ring elements.
Perm col
column permutation matrices
void solveUpdateRight(Real *vec)
int * orig
orig[p] original index from p
int solveRight4update(Real *vec, int *nonz, Real eps, Real *rhs, Real *forest, int *forestNum, int *forestIdx)
void vSolveLright(Real *vec, int *ridx, int &rn, Real eps)
Real lMemMult
factor of minimum Memory * number of nonzeros
void eliminateRowSingletons()
virtual void start()=0
start timer, resume accounting user, system and real time.
Real * work
Working array: must always be left as 0!
void forestUpdate(int col, Real *work, int num, int *nonz)
Performs the Forrest-Tomlin update of the LU factorization.
virtual Real stop()=0
stop timer, return accounted user time.
void spx_alloc(T &p, int n=1)
Allocate memory.
void solveLleftForestNoNZ(Real *vec)
static const Real verySparseFactor4right
int * len
used nonzeros per column vector
#define SOPLEX_FACTOR_MARKER
void solveUpdateLeft2(Real *vec1, Real *vec2)
int & index(int n)
Reference to index of n 'th nonzero.
int solveLleftEps(Real *vec, int *nonz, Real eps)
int size
size of arrays val and idx
void solveUright2(Real *work1, Real *vec1, Real *work2, Real *vec2)
SLinSolver::Status stat
Status indicator.
Pring * pivot_row
row index handlers for Real linked list
#define MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
int used
used entries of arrays idx and val
void solveLright(Real *vec)
void vSolveUpdateRightNoNZ(Real *vec, Real)
int startSize
size of array start
void vSolveLright2(Real *vec, int *ridx, int &rn, Real eps, Real *vec2, int *ridx2, int &rn2, Real eps2)
Pring * pivot_colNZ
lists for columns to number of nonzeros
static Real epsilonUpdate()
int factorCount
Number of factorizations.
Real * val
hold nonzero values: this is only initialized in the end of the factorization with DEFAULT updates...
int stage
stage of the structure
int pos
position of pivot column in row
int firstUnused
number of first unused L vector
Temp temp
Temporary storage.
void setPivot(const int p_stage, const int p_col, const int p_row, const Real val)
void solveUpdateLeft(Real *vec)
void eliminateNucleus(const Real eps, const Real threshold)
void vSolveRight4update2sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real eps2, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2, Real *forest, int *forestNum, int *forestIdx)
sparse version of above method
int thedim
dimension of factorized matrix
struct soplex::CLUFactor::U::Row row
int * start
starting positions in val and idx
Debugging, floating point type and parameter definitions.
int size
size of arrays val and idx
void spx_realloc(T &p, int n)
Change amount of allocated memory.
int * idx
indices of L vectors
void factor(const SVector **vec, Real threshold, Real eps)
epsilon for zero detection
int * perm
perm[i] permuted index from i
bool isConsistent() const
void eliminateColSingletons()
void solveUleft2(Real *work1, Real *vec1, Real *work2, Real *vec2)
Everything should be within this namespace.
void vSolveRightNoNZ(Real *vec, Real eps, Real *rhs, int *ridx, int rn)
void remaxRow(int p_row, int len)
int vSolveLeft(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn)
int solveLleft2forest(Real *vec1, int *, Real *vec2, Real)
void vSolveLeftNoNZ(Real eps, Real *vec, Real *rhs, int *ridx, int rn)
int firstUpdate
number of first update L vector
Pring pivots
ring of selected pivot rows
int vSolveLeft3(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real *rhs2, int *ridx2, int rn2, Real *vec3, Real *rhs3, int *ridx3, int rn3)
Real rowMemMult
factor of minimum Memory * number of nonzeros
void solveLright2(Real *vec1, Real *vec2)
Dring list
Double linked ringlist of vector indices in the order they appear in the row file.
int vSolveRight4update2(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int rn, Real *vec2, Real eps2, Real *rhs2, int *ridx2, int rn2, Real *forest, int *forestNum, int *forestIdx)
Timer * factorTime
Time spent in factorizations.
void clear()
clears the structure
int solveLeft2(Real *vec1, int *nonz, Real *vec2, Real eps, Real *rhs1, Real *rhs2)
void init(int p_dim)
initialization
static int deQueueMax(int *heap, int *size)
int solveLeftEps(Real *vec, Real *rhs, int *nonz, Real eps)
void vSolveUrightNoNZ(Real *vec, Real *rhs, int *ridx, int rn, Real eps)
void forestReMaxCol(int col, int len)
Implementation of sparse LU factorization.
int * idx
hold row indices of nonzeros
int idx
index of pivot row
void eliminatePivot(int prow, int pos, Real eps)
void remaxCol(int p_col, int len)
int * ridx
indices of rows of L
int vSolveUright2(Real *vec, int *vidx, Real *rhs, int *ridx, int rn, Real eps, Real *vec2, Real *rhs2, int *ridx2, int rn2, Real eps2)
Real colMemMult
factor of minimum Memory * number of nonzeros
int updateRow(int r, int lv, int prow, int pcol, Real pval, Real eps)
int vSolveUpdateRight(Real *vec, int *ridx, int n, Real eps)
int * max
maximum available nonzeros per colunn: start[i] + max[i] == start[elem[i].next->idx] len[i] <= max[i]...
static const Real verySparseFactor
int makeLvec(int p_len, int p_row)
void updateSolutionVectorLright(Real change, int j, Real &vec, int *idx, int &nnz)
int * row
column indices of L vectors
void vSolveLeft2sparse(Real eps, Real *vec, int *idx, Real *rhs, int *ridx, int &rn, Real *vec2, int *idx2, Real *rhs2, int *ridx2, int &rn2)
sparse version of solving 2 systems of equations
static Real epsilonPivot()
int * rbeg
start of rows in rval and ridx
int * s_cact
lengths of columns of active submatrix
void solveRight(Real *vec, Real *rhs)
void initFactorMatrix(const SVector **vec, const Real eps)
void solveLleft2(Real *vec1, int *, Real *vec2, Real)
static void enQueueMax(int *heap, int *size, int elem)
Exception class for status exceptions during the computationsThis class is derived from the SoPlex ex...
The loaded matrix is singular.
int nzCnt
number of nonzeros in U
void selectPivots(Real threshold)
Real * val
hold nonzero values
bool isZero(Real a, Real eps=Param::epsilon())
returns true iff |a| <= eps
int * rperm
original row permutation
int * start
starting positions in val and idx
int * idx
hold column indices of nonzeros
int used
used entries of array idx
void solveRight2(Real *vec1, Real *vec2, Real *rhs1, Real *rhs2)
int vSolveUright(Real *vec, int *vidx, Real *rhs, int *ridx, int rn, Real eps)
int solveUrightEps(Real *vec, int *nonz, Real eps, Real *rhs)
Real * val
values of L vectors
void spx_free(T &p)
Release memory.
void solveLleftNoNZ(Real *vec)
static int deQueueMin(int *heap, int *size)