|
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] )
237 if ( pivot_col != 0 )
240 if ( pivot_colNZ != 0 )
243 if ( pivot_row != 0 )
246 if ( pivot_rowNZ != 0 )
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 )
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++ );
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;
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++ );
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 );
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 );
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 );
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 )
|