|
Go to the documentation of this file. 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] ) 134 if ( elem < heap[i] ) 147 for ( i = 1; i < *size; ++i ) 148 for ( j = 0; j < i; ++j ) 149 assert( heap[i] != heap[j] ); 161 e = heap[s = --( *size )]; 164 for ( j = 0, i = 1; i < s; i = 2 * j + 1 ) 197 if ( i < *size && e > heap[i] ) 259 for ( int i = 0; i < thedim; ++i ) 277 diag[p_row] = 1.0 / val; 301 for ( ring = list-> next; ring != list; ring = ring-> next ) 305 if ( l_rbeg[l_row] != n ) 311 assert( l_rlen[l_row] <= l_rmax[l_row] ); 313 l_rmax[l_row] = l_rlen[l_row]; 314 j = i + l_rlen[l_row]; 316 for ( ; i < j; ++i, ++n ) 319 l_ridx[n] = l_ridx[i]; 320 l_rval[n] = l_rval[i]; 325 while ( ring != list ); 327 goto terminatePackRows; 332 l_rmax[l_row] = l_rlen[l_row]; 359 for ( ring = list-> next; ring != list; ring = ring-> next ) 363 if ( cbeg[colno] != n ) 370 cmax[colno] = clen[colno]; 381 while ( ring != list ); 383 goto terminatePackColumns; 388 cmax[colno] = clen[colno]; 391 terminatePackColumns : 402 assert( u. row. max[p_row] < len ); 406 int delta = len - u. row. max[p_row]; 411 delta = len - u. row. max[p_row]; 421 && "ERROR: could not allocate memory for row file" ); 444 && "ERROR: could not allocate memory for row file" ); 459 for ( ; i < k; ++i, ++j ) 488 for ( ring = list-> next; ring != list; ring = ring-> next ) 492 if ( l_cbeg[l_col] != n ) 499 l_cmax[l_col] = l_clen[l_col]; 500 j = i + l_clen[l_col]; 503 l_cidx[n++] = l_cidx[i]; 507 while ( ring != list ); 509 goto terminatePackColumns; 514 l_cmax[l_col] = l_clen[l_col]; 517 terminatePackColumns : 528 assert( u. col. max[p_col] < len ); 532 int delta = len - u. col. max[p_col]; 537 delta = len - u. col. max[p_col]; 547 && "ERROR: could not allocate memory for column file" ); 571 && "ERROR: could not allocate memory for column file" ); 597 assert( u. col. max[p_col] < len ); 601 int delta = len - u. col. max[p_col]; 606 delta = len - u. col. max[p_col]; 614 && "ERROR: could not allocate memory for column file" ); 637 && "ERROR: could not allocate memory for column file" ); 693 int i, j, k, h, m, n; 725 for ( i += j - 1; i >= j; --i ) 729 h = --( rlen[m] ) + k; 731 while ( ridx[k] != p_col ) 756 if ( num > cmax[p_col] ) 767 for ( j = 0; j < num; ++j ) 775 if ( spxAbs( x ) > l_maxabs ) 779 assert( k - cbeg[p_col] < cmax[p_col] ); 786 if ( rmax[i] <= rlen[i] ) 793 h = rbeg[i] + ( rlen[i] )++; 805 nzCnt += ( clen[p_col] = k - cbeg[p_col] ); 822 for ( i = 0; i < dim; ++i ) 829 if ( spxAbs( x ) > l_maxabs ) 835 clen[p_col] = k - cbeg[p_col]; 844 assert( k - cbeg[p_col] < cmax[p_col] ); 851 if ( rmax[i] <= rlen[i] ) 858 h = rbeg[i] + ( rlen[i] )++; 870 nzCnt += ( clen[p_col] = k - cbeg[p_col] ); 872 if ( cbeg[p_col] + cmax[p_col] == u. col. used ) 875 cmax[p_col] = clen[p_col]; 891 memmove(&rorig[c], &rorig[c + 1], ( unsigned int)(r - c) * sizeof( int)); 895 for ( i = c; i <= r; ++i ) 903 memmove(&corig[c], &corig[c + 1], ( unsigned int)(r - c) * sizeof( int)); 907 for ( i = c; i <= r; ++i ) 919 if ( i < verySparseFactor * ( dim - c ) ) 936 for ( i += j - 1; i >= j; --i ) 941 m = --( clen[k] ) + cbeg[k]; 943 for ( h = m; cidx[h] != rowno; --h ) 959 assert(( num == 0 ) || ( nonz != 0 ) ); 967 for ( i = 0; i < num; ++i ) 968 assert( p_work[corig[nonz[i]]] != 0 ); 978 assert( p_work[k] != 0 ); 982 x = p_work[k] * diag[n]; 992 if ( spxAbs( x ) > l_maxabs ) 1031 diag[rowno] = 1 / x; 1038 if ( rmax[rowno] < num ) 1052 for ( i = 0; i < num; ++i ) 1065 if ( spxAbs( x ) > l_maxabs ) 1076 if ( clen[j] >= cmax[j] ) 1083 cval[cbeg[j] + clen[j]] = x; 1085 cidx[cbeg[j] + clen[j]++] = rowno; 1089 rlen[rowno] = n - rbeg[rowno]; 1095 for ( i += j - 1; i >= j; --i ) 1098 p_work[k] = rval[i]; 1099 m = --( clen[k] ) + cbeg[k]; 1101 for ( h = m; cidx[h] != rowno; --h ) 1117 for ( i = c; i < r; ++i ) 1121 if ( p_work[k] != 0 ) 1124 x = p_work[k] * diag[n]; 1130 if ( spxAbs( x ) > l_maxabs ) 1137 for ( ; j < m; ++j ) 1138 p_work[ridx[j]] -= x * rval[j]; 1161 diag[rowno] = 1 / x; 1171 for ( i = r + 1; i < dim; ++i ) 1172 if ( p_work[corig[i]] != 0 ) 1175 if ( rmax[rowno] < n ) 1189 for ( i = r + 1; i < dim; ++i ) 1196 if ( spxAbs( x ) > l_maxabs ) 1207 if ( clen[j] >= cmax[j] ) 1214 cval[cbeg[j] + clen[j]] = x; 1216 cidx[cbeg[j] + clen[j]++] = rowno; 1220 rlen[rowno] = n - rbeg[rowno]; 1231 i = rbeg[rowno] + --( rlen[rowno] ); 1232 diag[rowno] = 1 / rval[i]; 1234 for ( j = i = --( clen[p_col] ) + cbeg[p_col]; cidx[i] != rowno; --i ) 1259 assert( p_work[p_col] != 0 ); 1260 rezi = 1 / p_work[p_col]; 1268 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i ) 1271 lval[ll] = rezi * p_work[j]; 1278 lval[ll] = 1 - rezi; 1281 for ( --i; i >= 0; --i ) 1285 lval[ll] = x = rezi * p_work[j]; 1305 assert( p_work[p_col] != 0 ); 1306 rezi = 1 / p_work[p_col]; 1312 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i ) 1315 lval[ll] = rezi * p_work[j]; 1321 lval[ll] = 1 - rezi; 1324 for ( --i; i >= 0; --i ) 1328 lval[ll] = x = rezi * p_work[j]; 1370 Dring *rring, *lastrring; 1371 Dring *cring, *lastcring; 1381 for ( int i = 0; i < thedim; i++ ) 1386 for ( int i = 0; i < thedim; i++ ) 1397 for ( int j = 0; j < k; ++j ) 1430 lastrring-> next = rring; 1438 lastcring-> next = cring; 1442 for ( int i = 0; i < thedim; i++ ) 1448 rring-> prev = lastrring; 1449 lastrring-> next = rring; 1454 cring-> prev = lastcring; 1455 lastcring-> next = cring; 1479 for ( int i = 0; i < thedim; i++ ) 1489 for ( int j = 0; j < psv-> size() && nnonzeros <= 1; j++ ) 1491 if ( psv-> value(j) != 0 ) 1496 if ( nnonzeros == 0 ) 1504 if ( nnonzeros == 1 ) 1510 for ( j = 0; psv-> value(j) == 0; j++ ) 1513 assert( j < psv->size() ); 1523 x = psv-> value( j ); 1546 assert( nnonzeros >= 2 ); 1549 for ( int j = 0; j < psv-> size(); j++ ) 1551 x = psv-> value( j ); 1556 k = psv-> index( j ); 1575 assert( nnonzeros >= 2 ); 1596 int p_col, p_row, newrow; 1612 assert( p_row >= 0 ); 1616 for ( j = 0; j < len; ++j ) 1625 for ( k = n; u. col. idx[k] != p_row; ++k ) 1643 if ( rperm[newrow] >= 0 ) 1653 for ( k = n; u. row. idx[k] != p_col; --k ) 1689 int p_row, p_col, len, rs, lk; 1698 for ( i = 0; i < thedim; ++i ) 1700 if ( rperm[i] < 0 && u. row. len[i] == 1 ) 1716 setPivot( rs, p_col, p_row, pval ); 1726 i = ( u. col. len[p_col] -= i ); 1728 for ( ; i < len; ++i ) 1739 for ( j = k; u. row. idx[j] != p_col; --j ) 1799 for ( i = 0; i < thedim; ++i ) 1884 for ( ; ( r = idx[i] ) != prow; ++i ) 1891 for ( j = k; u. row. idx[j] != pcol; --j ) 1922 assert( i < len && "ERROR: pivot column does not contain pivot row" ); 1924 for ( ++i; i < len; ++i ) 1932 for ( j = k; u. row. idx[j] != pcol; --j ) 1994 for ( i = j; ( c = u. row. idx[i] ) != pcol; --i ) 1998 for ( k = m; u. col. idx[k] != prow; ++k ) 2024 for ( --i; i >= j; --i ) 2029 for ( k = m; u. col. idx[k] != prow; ++k ) 2074 if ( candidates > 4 ) 2089 len = u. row. len[rw] + beg - 1; 2099 for ( i = len - 1; i >= beg; --i ) 2108 l_maxabs *= threshold; 2114 for ( i = len; i >= beg; --i ) 2120 if ( j < mkwtz && spxAbs( x ) > l_maxabs ) 2136 len = u. col. len[cl] + beg - 1; 2144 for ( i = len; i >= beg; --i ) 2177 for ( kk = m + u. row. len[k] - 1; kk >= m; --kk ) 2190 for ( --kk; kk > m; --kk ) 2199 l_maxabs *= threshold; 2201 if ( spxAbs( x ) > l_maxabs ) 2207 if ( j <= count + 1 ) 2236 if ( pr-> idx == rw || pr-> mkwtz >= mkwtz ) 2242 if ( pr-> idx != rw ) 2250 if ( num >= candidates ) 2282 int c, i, j, k, ll, m, n; 2291 for ( j = m; u. row. idx[j] != pcol; --j ) 2312 for ( j = m - 1; j >= n; --j ) 2341 for ( i = k; u. col. idx[i] != r; --i ) 2354 if ( ll + fill > u. row. max[r] ) 2411 int i, j, k, m = -1; 2416 int plen = --( u. row. len[prow] ); 2417 int pend = pbeg + plen; 2446 for ( i = pbeg; i < pend; ++i ) 2454 for ( k = m; u. col. idx[k] != prow; ++k ) 2479 for ( ++i; i < m; ++i ) 2491 for ( i = u. row. start[prow], pend = i + plen; i < pend; ++i ) 2525 for ( i = 0; i < thedim; ++i ) 2554 "ERROR: no pivot element selected" ); 2557 pivot = pivot-> next ) 2580 "ERROR: one row must be left" ); 2582 "ERROR: one col must be left" ); 2600 for ( i = 0; i < thedim; i++ ) 2605 for ( i = 0; i < thedim; i++ ) 2616 assert(( *idx >= 0 ) && ( *idx < thedim ) ); 2620 assert(( k >= 0 ) && ( k < u. col. size ) ); 2702 for ( i = thedim; i--; *l_rbeg++ = 0 ) 2704 *rorig++ = *rrorig++; 2705 *rperm++ = *rrperm++; 2710 l_rbeg = l. rbeg + 1; 2712 for ( i = mem; i--; ) 2717 for ( m = 0, i = thedim; i--; l_rbeg++ ) 2726 l_rbeg = l. rbeg + 1; 2728 for ( i = j = 0; i < vecs; ++i ) 2733 for ( ; j < beg[i + 1]; j++ ) 2735 k = l_rbeg[*idx++]++; 2738 l_rval[k] = val[validx]; 2745 assert( l. rbeg[0] == 0 ); 2757 MSG_DEBUG( std::cout << "CLUFactorRational::factor()\n" ); 2829 for ( i = 0; i < thedim; ++i ) 2832 std::cout << "DCLUFA01 diag[" << i << "]: [" << col. orig[ row. perm[i]] 2833 << "] = " << diag[i] << std::endl; 2835 for ( j = 0; j < u. row. len[i]; ++j ) 2836 std::cout << "DCLUFA02 u[" << i << "]: [" 2843 for ( i = 0; i < thedim; ++i ) 2848 std::cout << "DCLUFA03 l[" << i << "]" << std::endl; 2851 std::cout << "DCLUFA04 l[" << k - l. start[j] 2852 << "]: [" << l. idx[k] 2853 << "] = " << l. val[k] << std::endl; 2906 int newsize = int( 0.2 * l. val. dim() + size ); 2922 int* p_lrow = l. row; 2927 assert( p_len > 0 && "ERROR: no empty columns allowed in L vectors" ); 2943 #ifdef ENABLE_CONSISTENCY_CHECKS 2960 assert( ring-> idx >= 0 ); 2962 assert( ring-> prev-> next == ring ); 2985 assert( ring-> idx >= 0 ); 2987 assert( ring-> prev-> next == ring ); 3002 for ( i = 0; i < thedim; ++i ) 3011 for ( i = 0; i < thedim; ++i ) 3038 for ( i = 0; i < thedim; ++i ) 3060 for ( i = 0; i < thedim - temp. stage; ++i ) 3074 #endif // CONSISTENCY_CHECKS 3082 for ( int i = thedim - 1; i >= 0; i-- ) 3106 int *cidx, *clen, *cbeg; 3122 for ( i = thedim - 1; i >= 0; --i ) 3125 x = diag[r] * rhs[r]; 3132 val = &cval[cbeg[c]]; 3133 idx = &cidx[cbeg[c]]; 3137 rhs[*idx++] -= x * ( *val++ ); 3148 int *cidx, *clen, *cbeg; 3162 for ( i = thedim - 1; i >= 0; --i ) 3166 p_work1[c] = x1 = diag[r] * vec1[r]; 3167 p_work2[c] = x2 = diag[r] * vec2[r]; 3168 vec1[r] = vec2[r] = 0; 3170 if ( x1 != 0 && x2 != 0 ) 3172 val = &cval[cbeg[c]]; 3173 idx = &cidx[cbeg[c]]; 3178 vec1[*idx] -= x1 * ( *val ); 3179 vec2[*idx++] -= x2 * ( *val++ ); 3185 val = &cval[cbeg[c]]; 3186 idx = &cidx[cbeg[c]]; 3190 vec1[*idx++] -= x1 * ( *val++ ); 3195 val = &cval[cbeg[c]]; 3196 idx = &cidx[cbeg[c]]; 3200 vec2[*idx++] -= x2 * ( *val++ ); 3209 int *cidx, *clen, *cbeg; 3210 bool notzero1, notzero2; 3226 for ( i = thedim - 1; i >= 0; --i ) 3230 p_work1[c] = x1 = diag[r] * vec1[r]; 3231 p_work2[c] = x2 = diag[r] * vec2[r]; 3232 vec1[r] = vec2[r] = 0; 3233 notzero1 = x1 != 0 ? 1 : 0; 3234 notzero2 = x2 != 0 ? 1 : 0; 3236 if ( notzero1 && notzero2 ) 3240 val = &cval[cbeg[c]]; 3241 idx = &cidx[cbeg[c]]; 3246 vec1[*idx] -= x1 * ( *val ); 3247 vec2[*idx++] -= x2 * ( *val++ ); 3256 val = &cval[cbeg[c]]; 3257 idx = &cidx[cbeg[c]]; 3261 vec1[*idx++] -= x1 * ( *val++ ); 3267 val = &cval[cbeg[c]]; 3268 idx = &cidx[cbeg[c]]; 3272 vec2[*idx++] -= x2 * ( *val++ ); 3290 int *lrow, *lidx, *idx; 3300 for ( i = 0; i < end; ++i ) 3302 if (( x = vec[lrow[i]] ) != 0 ) 3307 MSG_DEBUG( std::cout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl; ) 3313 for ( j = lbeg[i + 1]; j > k; --j ) 3315 MSG_DEBUG( std::cout << " -> y" << *idx << " -= " << x << " * " << *val << " = " << x * ( *val ) << " -> " << vec[*idx] - x * ( *val ) << std::endl; ) 3316 vec[*idx++] -= x * ( *val++ ); 3323 MSG_DEBUG( std::cout << "performing FT updates..." << std::endl; ) 3327 for ( ; i < end; ++i ) 3334 for ( j = lbeg[i + 1]; j > k; --j ) 3335 x += vec[*idx++] * ( *val++ ); 3339 MSG_DEBUG( std::cout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl; ) 3342 MSG_DEBUG( std::cout << "finished FT updates." << std::endl; ) 3353 int *lrow, *lidx, *idx; 3363 for ( i = 0; i < end; ++i ) 3368 if ( x1 != 0 && x2 != 0 ) 3374 for ( j = lbeg[i + 1]; j > k; --j ) 3376 vec1[*idx] -= x1 * ( *val ); 3377 vec2[*idx++] -= x2 * ( *val++ ); 3387 for ( j = lbeg[i + 1]; j > k; --j ) 3388 vec1[*idx++] -= x1 * ( *val++ ); 3397 for ( j = lbeg[i + 1]; j > k; --j ) 3398 vec2[*idx++] -= x2 * ( *val++ ); 3406 for ( ; i < end; ++i ) 3414 for ( j = lbeg[i + 1]; j > k; --j ) 3416 x1 += vec1[*idx] * ( *val ); 3417 x2 += vec2[*idx++] * ( *val++ ); 3420 vec1[lrow[i]] -= x1; 3422 vec2[lrow[i]] -= x2; 3433 int *lrow, *lidx, *idx; 3447 if (( x = vec[lrow[i]] ) != 0 ) 3453 for ( j = lbeg[i + 1]; j > k; --j ) 3454 vec[*idx++] -= x * ( *val++ ); 3484 if ( x1 != 0 && x2 != 0 ) 3490 for ( j = lbeg[i + 1]; j > k; --j ) 3492 vec1[*idx] -= x1 * ( *val ); 3493 vec2[*idx++] -= x2 * ( *val++ ); 3503 for ( j = lbeg[i + 1]; j > k; --j ) 3504 vec1[*idx++] -= x1 * ( *val++ ); 3513 for ( j = lbeg[i + 1]; j > k; --j ) 3514 vec2[*idx++] -= x2 * ( *val++ ); 3528 for ( int i = 0; i < thedim; i++ ) 3532 n += rhs[i] != 0 ? 1 : 0; 3572 for ( int i = 0; i < thedim; i++ ) 3575 forest[i] = rhs1[i]; 3576 n += rhs1[i] != 0 ? 1 : 0; 3612 for ( int i = 0; i < thedim; ++i ) 3625 #if defined(WITH_WARNINGS) || !defined(NDEBUG) 3645 for ( int m = u. row. start[r]; m < end; m++ ) 3666 int *ridx, *rlen, *rbeg, *idx; 3677 for ( i = 0; i < thedim; ++i ) 3684 if (( x1 != 0 ) && ( x2 != 0 ) ) 3694 for ( int m = rlen[r]; m != 0; --m ) 3696 vec1[*idx] -= x1 * ( *val ); 3697 vec2[*idx] -= x2 * ( *val++ ); 3710 for ( int m = rlen[r]; m != 0; --m ) 3711 vec1[*idx++] -= x1 * ( *val++ ); 3722 for ( int m = rlen[r]; m != 0; --m ) 3723 vec2[*idx++] -= x2 * ( *val++ ); 3736 int *lidx, *idx, *lrow; 3760 for ( j = lbeg[i + 1]; j > k; --j ) 3762 vec1[*idx] -= x1 * ( *val ); 3763 vec2[*idx++] -= x2 * ( *val++ ); 3772 for ( j = lbeg[i + 1]; j > k; --j ) 3773 vec1[*idx++] -= x1 * ( *val++ ); 3783 for ( j = lbeg[i + 1]; j > k; --j ) 3784 vec2[*idx++] -= x2 * ( *val++ ); 3815 for ( ; i >= 0; --i ) 3823 for ( j = lbeg[i + 1]; j > k; --j ) 3825 x1 += vec1[*idx] * ( *val ); 3826 x2 += vec2[*idx++] * ( *val++ ); 3829 vec1[lrow[i]] -= x1; 3831 vec2[lrow[i]] -= x2; 3840 x1not0 = ( x1 != 0 ); 3841 x2not0 = ( x2 != 0 ); 3843 if ( x1not0 && x2not0 ) 3846 j = rbeg[r + 1] - k; 3853 vec1[*idx] -= x1 * *val; 3854 vec2[*idx++] -= x2 * *val++; 3861 j = rbeg[r + 1] - k; 3868 vec1[*idx++] -= x1 * *val++; 3875 j = rbeg[r + 1] - k; 3882 vec2[*idx++] -= x2 * *val++; 3895 int *idx, *lidx, *lrow, *lbeg; 3906 if (( x = vec[lrow[i]] ) != 0 ) 3915 for ( j = lbeg[i + 1]; j > k; --j ) 3916 vec[*idx++] -= x * ( *val++ ); 3944 for ( int j = lbeg[i + 1]; j > k; --j ) 3945 x += vec[*idx++] * ( *val++ ); 3951 for ( int i = thedim - 1; i >= 0; --i ) 3961 for ( int k = l. rbeg[r]; k < l. rbeg[r + 1]; k++ ) 3965 assert( l. rperm[j] < i ); 3967 vec[j] -= x * l. rval[k]; 3972 #endif // WITH_L_ROWS 4003 for ( j = lbeg[i + 1]; j > k; --j ) 4004 x += vec[*idx++] * ( *val++ ); 4020 j = rbeg[r + 1] - k; 4027 vec[*idx++] -= x * *val++; 4044 int *lrow, *lidx, *idx; 4063 for ( j = lbeg[i + 1]; j > k; --j ) 4064 x += vec[*idx++] * ( *val++ ); 4075 int *lrow, *lidx, *idx; 4095 for ( j = lbeg[i + 1]; j > k; --j ) 4097 x1 += vec1[*idx] * ( *val ); 4098 x2 += vec2[*idx++] * ( *val++ ); 4101 vec1[lrow[i]] -= x1; 4103 vec2[lrow[i]] -= x2; 4112 int *lrow, *lidx, *idx; 4127 assert( k >= 0 && k < l. val. dim() ); 4132 for ( j = lbeg[i + 1]; j > k; --j ) 4134 assert( *idx >= 0 && *idx < thedim ); 4135 x += vec[*idx++] * ( *val++ ); 4222 Rational* rhs, int* rhsidx, int rhsn ) 4225 int i, j, k, n, r, c; 4226 int *rorig, *corig, *cperm; 4227 int *ridx, *rlen, *rbeg, *idx; 4237 for ( i = 0; i < rhsn; ) 4253 assert( i >= 0 && i < thedim ); 4255 assert( c >= 0 && c < thedim ); 4262 assert( r >= 0 && r < thedim ); 4271 for ( int m = rlen[r]; m; --m ) 4274 assert( j >= 0 && j < thedim ); 4279 y = -x * ( *val++ ); 4289 y -= x * ( *val++ ); 4305 int *rorig, *corig, *cperm; 4306 int *ridx, *rlen, *rbeg, *idx; 4316 for ( i = 0; i < rhsn; ) 4330 assert( i >= 0 && i < thedim ); 4332 assert( c >= 0 && c < thedim ); 4339 assert( r >= 0 && r < thedim ); 4347 for ( int m = rlen[r]; m; --m ) 4350 assert( j >= 0 && j < thedim ); 4355 y = -x * ( *val++ ); 4365 y -= x * ( *val++ ); 4380 int *idx, *lidx, *lrow, *lbeg; 4390 assert( i >= 0 && i < l. val. dim() ); 4392 if (( x = vec[lrow[i]] ) != 0 ) 4395 assert( k >= 0 && k < l. val. dim() ); 4399 for ( j = lbeg[i + 1]; j > k; --j ) 4402 assert( m >= 0 && m < thedim ); 4407 y = -x * ( *val++ ); 4417 y -= x * ( *val++ ); 4433 int *idx, *lidx, *lrow, *lbeg; 4443 if (( x = vec[lrow[i]] ) != 0 ) 4445 assert( i >= 0 && i < l. val. dim() ); 4447 assert( k >= 0 && k < l. val. dim() ); 4451 for ( j = lbeg[i + 1]; j > k; --j ) 4453 assert( *idx >= 0 && *idx < thedim ); 4454 vec[*idx++] -= x * ( *val++ ); 4481 #pragma warn "Not yet implemented, define WITH_L_ROWS" 4487 for ( ; i >= 0; --i ) 4494 for ( j = lbeg[i + 1]; j > k; --j ) 4495 x += vec[*idx++] * ( *val++ ); 4503 for ( i = 0; i < rn; ) 4519 j = rbeg[r + 1] - k; 4525 assert( l. rperm[*idx] < i ); 4551 for ( i = 0; i < n; ++i ) 4584 for ( ; i >= 0; --i ) 4587 assert( k >= 0 && k < l. val. dim() ); 4592 for ( j = lbeg[i + 1]; j > k; --j ) 4594 assert( *idx >= 0 && *idx < thedim ); 4595 x += vec[*idx++] * ( *val++ ); 4610 j = rbeg[r + 1] - k; 4616 assert( l. rperm[*idx] < i ); 4617 vec[*idx++] -= x * *val++; 4631 int *lrow, *lidx, *idx; 4641 for ( i = 0; i < end; ++i ) 4651 for ( j = lbeg[i + 1]; j > k; --j ) 4653 assert( *idx >= 0 && *idx < thedim ); 4654 ridx[rn] = n = *idx++; 4655 rn += ( vec[n] == 0 ) ? 1 : 0; 4656 vec[n] -= x * ( *val++ ); 4666 for ( ; i < end; ++i ) 4673 for ( j = lbeg[i + 1]; j > k; --j ) 4675 assert( *idx >= 0 && *idx < thedim ); 4676 x += vec[*idx++] * ( *val++ ); 4679 ridx[rn] = j = lrow[i]; 4681 rn += ( vec[j] == 0 ) ? 1 : 0; 4691 Rational* vec2, int* ridx2, int* rn2ptr ) 4698 int *lrow, *lidx, *idx; 4711 for ( i = 0; i < end; ++i ) 4725 for ( j = lbeg[i + 1]; j > k; --j ) 4727 assert( *idx >= 0 && *idx < thedim ); 4728 ridx[rn] = ridx2[rn2] = n = *idx++; 4731 rn += ( y == 0 ) ? 1 : 0; 4732 rn2 += ( y2 == 0 ) ? 1 : 0; 4734 y2 -= x2 * ( *val++ ); 4743 for ( j = lbeg[i + 1]; j > k; --j ) 4745 assert( *idx >= 0 && *idx < thedim ); 4746 ridx[rn] = n = *idx++; 4748 rn += ( y == 0 ) ? 1 : 0; 4749 y -= x * ( *val++ ); 4762 for ( j = lbeg[i + 1]; j > k; --j ) 4764 assert( *idx >= 0 && *idx < thedim ); 4765 ridx2[rn2] = n = *idx++; 4767 rn2 += ( y2 == 0 ) ? 1 : 0; 4768 y2 -= x2 * ( *val++ ); 4779 for ( ; i < end; ++i ) 4786 for ( j = lbeg[i + 1]; j > k; --j ) 4788 assert( *idx >= 0 && *idx < thedim ); 4789 x += vec[*idx] * ( *val ); 4790 x2 += vec2[*idx++] * ( *val++ ); 4793 ridx[rn] = ridx2[rn2] = j = lrow[i]; 4795 rn += ( vec[j] == 0 ) ? 1 : 0; 4796 rn2 += ( vec2[j] == 0 ) ? 1 : 0; 4810 Rational* vec2, int* ridx2, int* rn2ptr, 4811 Rational* vec3, int* ridx3, int* rn3ptr) 4819 int *lrow, *lidx, *idx; 4833 for ( i = 0; i < end; ++i ) 4851 for ( j = lbeg[i + 1]; j > k; --j ) 4853 assert( *idx >= 0 && *idx < thedim ); 4854 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++; 4858 rn += ( y == 0 ) ? 1 : 0; 4859 rn2 += ( y2 == 0 ) ? 1 : 0; 4860 rn3 += ( y3 == 0 ) ? 1 : 0; 4862 y2 -= x2 * ( *val ); 4863 y3 -= x3 * ( *val++ ); 4875 for ( j = lbeg[i + 1]; j > k; --j ) 4877 assert( *idx >= 0 && *idx < thedim ); 4878 ridx[rn] = ridx2[rn2] = n = *idx++; 4881 rn += ( y == 0 ) ? 1 : 0; 4882 rn2 += ( y2 == 0 ) ? 1 : 0; 4884 y2 -= x2 * ( *val++ ); 4896 for ( j = lbeg[i + 1]; j > k; --j ) 4898 assert( *idx >= 0 && *idx < thedim ); 4899 ridx[rn] = ridx3[rn3] = n = *idx++; 4902 rn += ( y == 0 ) ? 1 : 0; 4903 rn3 += ( y3 == 0 ) ? 1 : 0; 4905 y3 -= x3 * ( *val++ ); 4915 for ( j = lbeg[i + 1]; j > k; --j ) 4917 assert( *idx >= 0 && *idx < thedim ); 4918 ridx[rn] = n = *idx++; 4920 rn += ( y == 0 ) ? 1 : 0; 4921 y -= x * ( *val++ ); 4937 for ( j = lbeg[i + 1]; j > k; --j ) 4939 assert( *idx >= 0 && *idx < thedim ); 4940 ridx2[rn2] = ridx3[rn3] = n = *idx++; 4943 rn2 += ( y2 == 0 ) ? 1 : 0; 4944 rn3 += ( y3 == 0 ) ? 1 : 0; 4945 y2 -= x2 * ( *val ); 4946 y3 -= x3 * ( *val++ ); 4956 for ( j = lbeg[i + 1]; j > k; --j ) 4958 assert( *idx >= 0 && *idx < thedim ); 4959 ridx2[rn2] = n = *idx++; 4961 rn2 += ( y2 == 0 ) ? 1 : 0; 4962 y2 -= x2 * ( *val++ ); 4976 for ( j = lbeg[i + 1]; j > k; --j ) 4978 assert( *idx >= 0 && *idx < thedim ); 4979 ridx3[rn3] = n = *idx++; 4981 rn3 += ( y3 == 0 ) ? 1 : 0; 4982 y3 -= x3 * ( *val++ ); 4993 for ( ; i < end; ++i ) 5000 for ( j = lbeg[i + 1]; j > k; --j ) 5002 assert( *idx >= 0 && *idx < thedim ); 5003 x += vec[*idx] * ( *val ); 5004 x2 += vec2[*idx] * ( *val ); 5005 x3 += vec3[*idx++] * ( *val++ ); 5008 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i]; 5010 rn += ( vec[j] == 0 ) ? 1 : 0; 5011 rn2 += ( vec2[j] == 0 ) ? 1 : 0; 5012 rn3 += ( vec3[j] == 0 ) ? 1 : 0; 5031 int i, j, k, r, c, n; 5034 int *cidx, *clen, *cbeg; 5056 assert( i >= 0 && i < thedim ); 5058 assert( r >= 0 && r < thedim ); 5060 x = diag[r] * rhs[r]; 5066 assert( c >= 0 && c < thedim ); 5069 val = &cval[cbeg[c]]; 5070 idx = &cidx[cbeg[c]]; 5075 assert( *idx >= 0 && *idx < thedim ); 5077 assert( k >= 0 && k < thedim ); 5082 y = -x * ( *val++ ); 5092 y -= x * ( *val++ ); 5098 if ( rn > i*verySparseFactor4right ) 5100 for ( i = *ridx; i >= 0; --i ) 5103 assert( r >= 0 && r < thedim ); 5104 x = diag[r] * rhs[r]; 5110 assert( c >= 0 && c < thedim ); 5113 val = &cval[cbeg[c]]; 5114 idx = &cidx[cbeg[c]]; 5119 assert( *idx >= 0 && *idx < thedim ); 5120 rhs[*idx++] -= x * ( *val++ ); 5139 int *cidx, *clen, *cbeg; 5156 if ( rn > *ridx * verySparseFactor4right ) 5158 for ( i = *ridx; i >= 0; --i ) 5160 assert( i >= 0 && i < thedim ); 5162 assert( r >= 0 && r < thedim ); 5163 x = diag[r] * rhs[r]; 5170 val = &cval[cbeg[c]]; 5171 idx = &cidx[cbeg[c]]; 5176 assert( *idx >= 0 && *idx < thedim ); 5177 rhs[*idx++] -= x * ( *val++ ); 5189 assert( i >= 0 && i < thedim ); 5193 assert( r >= 0 && r < thedim ); 5195 x = diag[r] * rhs[r]; 5203 val = &cval[cbeg[c]]; 5204 idx = &cidx[cbeg[c]]; 5210 assert( k >= 0 && k < thedim ); 5215 y = -x * ( *val++ ); 5225 y -= x * ( *val++ ); 5238 int i, j, k, r, c, n; 5241 int *cidx, *clen, *cbeg; 5259 while ( rn + rn2 > 0 ) 5269 if ( *ridx2 > *ridx ) 5272 if ( *ridx2 < *ridx ) 5280 assert( i >= 0 && i < thedim ); 5283 assert( r >= 0 && r < thedim ); 5285 x = diag[r] * rhs[r]; 5286 x2 = diag[r] * rhs2[r]; 5296 val = &cval[cbeg[c]]; 5297 idx = &cidx[cbeg[c]]; 5305 assert( k >= 0 && k < thedim ); 5310 y2 = -x2 * ( *val ); 5320 y2 -= x2 * ( *val ); 5329 y = -x * ( *val++ ); 5339 y -= x * ( *val++ ); 5350 assert( k >= 0 && k < thedim ); 5355 y = -x * ( *val++ ); 5365 y -= x * ( *val++ ); 5376 assert( c >= 0 && c < thedim ); 5378 val = &cval[cbeg[c]]; 5379 idx = &cidx[cbeg[c]]; 5385 assert( k >= 0 && k < thedim ); 5390 y2 = -x2 * ( *val++ ); 5400 y2 -= x2 * ( *val++ ); 5407 if ( rn + rn2 > i*verySparseFactor4right ) 5409 if ( *ridx > *ridx2 ) 5414 for ( ; i >= 0; --i ) 5418 assert( r >= 0 && r < thedim ); 5419 x = diag[r] * rhs[r]; 5420 x2 = diag[r] * rhs2[r]; 5427 assert( c >= 0 && c < thedim ); 5429 val = &cval[cbeg[c]]; 5430 idx = &cidx[cbeg[c]]; 5440 assert( *idx >= 0 && *idx < thedim ); 5441 rhs[*idx] -= x * ( *val ); 5442 rhs2[*idx++] -= x2 * ( *val++ ); 5449 assert( *idx >= 0 && *idx < thedim ); 5450 rhs2[*idx++] -= x2 * ( *val++ ); 5458 assert( c >= 0 && c < thedim ); 5461 val = &cval[cbeg[c]]; 5462 idx = &cidx[cbeg[c]]; 5467 assert( *idx >= 0 && *idx < thedim ); 5468 rhs[*idx++] -= x * ( *val++ ); 5486 int *lrow, *lidx, *idx; 5499 assert( i >= 0 && i < thedim ); 5505 assert( k >= 0 && k < l. val. dim() ); 5509 for ( j = lbeg[i + 1]; j > k; --j ) 5511 int m = ridx[n] = *idx++; 5512 assert( m >= 0 && m < thedim ); 5514 n += ( y == 0 ) ? 1 : 0; 5515 y = y - x * ( *val++ ); 5531 int *lrow, *lidx, *idx; 5544 assert( i >= 0 && i < thedim ); 5546 if (( x = vec[lrow[i]] ) != 0 ) 5549 assert( k >= 0 && k < l. val. dim() ); 5553 for ( j = lbeg[i + 1]; j > k; --j ) 5555 assert( *idx >= 0 && *idx < thedim ); 5556 vec[*idx++] -= x * ( *val++ ); 5565 Rational* forest, int* forestNum, int* forestIdx ) 5577 int* it = forestIdx; 5581 for ( i = j = 0; i < rn; ++i ) 5584 assert( k >= 0 && k < thedim ); 5596 *forestNum = rn = j; 5606 for ( i = j = 0; i < rn; ++i ) 5609 assert( k >= 0 && k < thedim ); 5632 Rational* rhs2, int* ridx2, int rn2, 5633 Rational* forest, int* forestNum, int* forestIdx ) 5645 int* it = forestIdx; 5649 for ( i = j = 0; i < rn; ++i ) 5652 assert( k >= 0 && k < thedim ); 5664 *forestNum = rn = j; 5674 for ( i = j = 0; i < rn; ++i ) 5677 assert( k >= 0 && k < thedim ); 5689 if ( rn2 > thedim*verySparseFactor4right ) 5704 for ( i = j = 0; i < rn2; ++i ) 5707 assert( k >= 0 && k < thedim ); 5743 Rational* rhs2, int* ridx2, int rn2, 5745 Rational* rhs3, int* ridx3, int rn3, 5746 Rational* forest, int* forestNum, int* forestIdx ) 5749 vSolveLright3( rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3 ); 5750 assert( rn >= 0 && rn <= thedim ); 5751 assert( rn2 >= 0 && rn2 <= thedim ); 5752 assert( rn3 >= 0 && rn3 <= thedim ); 5762 int* it = forestIdx; 5766 for ( i = j = 0; i < rn; ++i ) 5769 assert( k >= 0 && k < thedim ); 5781 *forestNum = rn = j; 5791 for ( i = j = 0; i < rn; ++i ) 5794 assert( k >= 0 && k < thedim ); 5806 if ( rn2 > thedim*verySparseFactor4right ) 5818 for ( i = j = 0; i < rn2; ++i ) 5821 assert( k >= 0 && k < thedim ); 5835 if ( rn3 > thedim*verySparseFactor4right ) 5847 for ( i = j = 0; i < rn3; ++i ) 5850 assert( k >= 0 && k < thedim ); 5880 Rational* rhs2, int* ridx2, int rn2 ) 5883 assert( rn2 >= 0 && rn2 <= thedim ); 5885 if ( rn2 > thedim*verySparseFactor4right ) 5899 for ( i = j = 0; i < rn2; ++i ) 5902 assert( k >= 0 && k < thedim ); 5950 Rational* rhs2, int* ridx2, int rn2 ) 5978 Rational* rhs2, int* ridx2, int rn2, 5980 Rational* rhs3, int* ridx3, int rn3 ) 6011 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!
static const Real verySparseFactor4left
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)
#define init2DR(elem, ring)
Implementation of sparse LU factorization with Rational precision.
Exception classes for SoPlex.
int thedim dimension of factorized matrix
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)
int size() const Number of used indices.
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)
int dim() const Dimension of vector.
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
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.
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 *)
bool isConsistent() const
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.
|