|
Go to the documentation of this file.
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] )
238 if ( pivot_col != 0 )
241 if ( pivot_colNZ != 0 )
244 if ( pivot_row != 0 )
247 if ( pivot_rowNZ != 0 )
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 )
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;
3930 for ( i = thedim; i--; )
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++ );
4096 for ( i = thedim; i--; )
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 )
4668 assert( i < thedim );
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++ );
4688 for ( i = thedim; i--; )
4696 j = rbeg[r + 1] - k;
4702 assert( l. rperm[*idx] < i );
4703 vec[*idx++] -= x * *val++;
4717 int *lrow, *lidx, *idx;
4727 for ( i = 0; i < end; ++i )
4737 for ( j = lbeg[i + 1]; j > k; --j )
4739 assert( *idx >= 0 && *idx < thedim );
4740 ridx[rn] = n = *idx++;
4741 rn += ( vec[n] == 0 ) ? 1 : 0;
4742 vec[n] -= x * ( *val++ );
4752 for ( ; i < end; ++i )
4759 for ( j = lbeg[i + 1]; j > k; --j )
4761 assert( *idx >= 0 && *idx < thedim );
4762 x += vec[*idx++] * ( *val++ );
4765 ridx[rn] = j = lrow[i];
4767 rn += ( vec[j] == 0 ) ? 1 : 0;
4777 Real* vec, int* ridx, int* rnptr, Real eps,
4778 Real* vec2, int* ridx2, int* rn2ptr, Real eps2 )
4785 int *lrow, *lidx, *idx;
4798 for ( i = 0; i < end; ++i )
4812 for ( j = lbeg[i + 1]; j > k; --j )
4814 assert( *idx >= 0 && *idx < thedim );
4815 ridx[rn] = ridx2[rn2] = n = *idx++;
4818 rn += ( y == 0 ) ? 1 : 0;
4819 rn2 += ( y2 == 0 ) ? 1 : 0;
4821 y2 -= x2 * ( *val++ );
4828 for ( j = lbeg[i + 1]; j > k; --j )
4830 assert( *idx >= 0 && *idx < thedim );
4831 ridx[rn] = n = *idx++;
4833 rn += ( y == 0 ) ? 1 : 0;
4834 y -= x * ( *val++ );
4846 for ( j = lbeg[i + 1]; j > k; --j )
4848 assert( *idx >= 0 && *idx < thedim );
4849 ridx2[rn2] = n = *idx++;
4851 rn2 += ( y2 == 0 ) ? 1 : 0;
4852 y2 -= x2 * ( *val++ );
4862 for ( ; i < end; ++i )
4869 for ( j = lbeg[i + 1]; j > k; --j )
4871 assert( *idx >= 0 && *idx < thedim );
4872 x += vec[*idx] * ( *val );
4873 x2 += vec2[*idx++] * ( *val++ );
4876 ridx[rn] = ridx2[rn2] = j = lrow[i];
4878 rn += ( vec[j] == 0 ) ? 1 : 0;
4879 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
4893 Real* vec, int* ridx, int* rnptr, Real eps,
4894 Real* vec2, int* ridx2, int* rn2ptr, Real eps2,
4895 Real* vec3, int* ridx3, int* rn3ptr, Real eps3 )
4903 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 );
4938 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4942 rn += ( y == 0 ) ? 1 : 0;
4943 rn2 += ( y2 == 0 ) ? 1 : 0;
4944 rn3 += ( y3 == 0 ) ? 1 : 0;
4946 y2 -= x2 * ( *val );
4947 y3 -= x3 * ( *val++ );
4956 for ( j = lbeg[i + 1]; j > k; --j )
4958 assert( *idx >= 0 && *idx < thedim );
4959 ridx[rn] = ridx2[rn2] = n = *idx++;
4962 rn += ( y == 0 ) ? 1 : 0;
4963 rn2 += ( y2 == 0 ) ? 1 : 0;
4965 y2 -= x2 * ( *val++ );
4975 for ( j = lbeg[i + 1]; j > k; --j )
4977 assert( *idx >= 0 && *idx < thedim );
4978 ridx[rn] = ridx3[rn3] = n = *idx++;
4981 rn += ( y == 0 ) ? 1 : 0;
4982 rn3 += ( y3 == 0 ) ? 1 : 0;
4984 y3 -= x3 * ( *val++ );
4992 for ( j = lbeg[i + 1]; j > k; --j )
4994 assert( *idx >= 0 && *idx < thedim );
4995 ridx[rn] = n = *idx++;
4997 rn += ( y == 0 ) ? 1 : 0;
4998 y -= x * ( *val++ );
5013 for ( j = lbeg[i + 1]; j > k; --j )
5015 assert( *idx >= 0 && *idx < thedim );
5016 ridx2[rn2] = ridx3[rn3] = n = *idx++;
5019 rn2 += ( y2 == 0 ) ? 1 : 0;
5020 rn3 += ( y3 == 0 ) ? 1 : 0;
5021 y2 -= x2 * ( *val );
5022 y3 -= x3 * ( *val++ );
5030 for ( j = lbeg[i + 1]; j > k; --j )
5032 assert( *idx >= 0 && *idx < thedim );
5033 ridx2[rn2] = n = *idx++;
5035 rn2 += ( y2 == 0 ) ? 1 : 0;
5036 y2 -= x2 * ( *val++ );
5049 for ( j = lbeg[i + 1]; j > k; --j )
5051 assert( *idx >= 0 && *idx < thedim );
5052 ridx3[rn3] = n = *idx++;
5054 rn3 += ( y3 == 0 ) ? 1 : 0;
5055 y3 -= x3 * ( *val++ );
5065 for ( ; i < end; ++i )
5072 for ( j = lbeg[i + 1]; j > k; --j )
5074 assert( *idx >= 0 && *idx < thedim );
5075 x += vec[*idx] * ( *val );
5076 x2 += vec2[*idx] * ( *val );
5077 x3 += vec3[*idx++] * ( *val++ );
5080 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
5082 rn += ( vec[j] == 0 ) ? 1 : 0;
5083 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
5084 rn3 += ( vec3[j] == 0 ) ? 1 : 0;
5101 Real* rhs, int* ridx, int rn, Real eps )
5103 int i, j, k, r, c, n;
5106 int *cidx, *clen, *cbeg;
5129 assert( i >= 0 && i < thedim );
5131 assert( r >= 0 && r < thedim );
5133 x = diag[r] * rhs[r];
5139 assert( c >= 0 && c < thedim );
5142 val = &cval[cbeg[c]];
5143 idx = &cidx[cbeg[c]];
5148 assert( *idx >= 0 && *idx < thedim );
5150 assert( k >= 0 && k < thedim );
5155 y = -x * ( *val++ );
5165 y -= x * ( *val++ );
5173 for ( i = *ridx; i >= 0; --i )
5176 assert( r >= 0 && r < thedim );
5177 x = diag[r] * rhs[r];
5183 assert( c >= 0 && c < thedim );
5186 val = &cval[cbeg[c]];
5187 idx = &cidx[cbeg[c]];
5192 assert( *idx >= 0 && *idx < thedim );
5193 rhs[*idx++] -= x * ( *val++ );
5208 Real* rhs, int* ridx, int rn, Real eps )
5213 int *cidx, *clen, *cbeg;
5233 for ( i = *ridx; i >= 0; --i )
5235 assert( i >= 0 && i < thedim );
5237 assert( r >= 0 && r < thedim );
5238 x = diag[r] * rhs[r];
5245 val = &cval[cbeg[c]];
5246 idx = &cidx[cbeg[c]];
5251 assert( *idx >= 0 && *idx < thedim );
5252 rhs[*idx++] -= x * ( *val++ );
5264 assert( i >= 0 && i < thedim );
5268 assert( r >= 0 && r < thedim );
5270 x = diag[r] * rhs[r];
5278 val = &cval[cbeg[c]];
5279 idx = &cidx[cbeg[c]];
5285 assert( k >= 0 && k < thedim );
5290 y = -x * ( *val++ );
5300 y -= x * ( *val++ );
5311 Real* vec, int* vidx, Real* rhs, int* ridx, int rn, Real eps,
5312 Real* vec2, Real* rhs2, int* ridx2, int rn2, Real eps2 )
5314 int i, j, k, r, c, n;
5317 int *cidx, *clen, *cbeg;
5336 while ( rn + rn2 > 0 )
5346 if ( *ridx2 > *ridx )
5349 if ( *ridx2 < *ridx )
5357 assert( i >= 0 && i < thedim );
5360 assert( r >= 0 && r < thedim );
5362 x = diag[r] * rhs[r];
5363 x2 = diag[r] * rhs2[r];
5373 val = &cval[cbeg[c]];
5374 idx = &cidx[cbeg[c]];
5382 assert( k >= 0 && k < thedim );
5387 y2 = -x2 * ( *val );
5397 y2 -= x2 * ( *val );
5405 y = -x * ( *val++ );
5415 y -= x * ( *val++ );
5426 assert( k >= 0 && k < thedim );
5431 y = -x * ( *val++ );
5441 y -= x * ( *val++ );
5452 assert( c >= 0 && c < thedim );
5454 val = &cval[cbeg[c]];
5455 idx = &cidx[cbeg[c]];
5461 assert( k >= 0 && k < thedim );
5466 y2 = -x2 * ( *val++ );
5476 y2 -= x2 * ( *val++ );
5484 if ( *ridx > *ridx2 )
5489 for ( ; i >= 0; --i )
5491 assert( i < thedim );
5493 assert( r >= 0 && r < thedim );
5494 x = diag[r] * rhs[r];
5495 x2 = diag[r] * rhs2[r];
5502 assert( c >= 0 && c < thedim );
5504 val = &cval[cbeg[c]];
5505 idx = &cidx[cbeg[c]];
5515 assert( *idx >= 0 && *idx < thedim );
5516 rhs[*idx] -= x * ( *val );
5517 rhs2[*idx++] -= x2 * ( *val++ );
5524 assert( *idx >= 0 && *idx < thedim );
5525 rhs2[*idx++] -= x2 * ( *val++ );
5533 assert( c >= 0 && c < thedim );
5536 val = &cval[cbeg[c]];
5537 idx = &cidx[cbeg[c]];
5542 assert( *idx >= 0 && *idx < thedim );
5543 rhs[*idx++] -= x * ( *val++ );
5561 int *lrow, *lidx, *idx;
5574 assert( i >= 0 && i < thedim );
5580 assert( k >= 0 && k < l. size );
5584 for ( j = lbeg[i + 1]; j > k; --j )
5586 int m = ridx[n] = *idx++;
5587 assert( m >= 0 && m < thedim );
5589 n += ( y == 0 ) ? 1 : 0;
5590 y = y - x * ( *val++ );
5605 int *lrow, *lidx, *idx;
5618 assert( i >= 0 && i < thedim );
5620 if (( x = vec[lrow[i]] ) != 0.0 )
5623 assert( k >= 0 && k < l. size );
5627 for ( j = lbeg[i + 1]; j > k; --j )
5629 assert( *idx >= 0 && *idx < thedim );
5630 vec[*idx++] -= x * ( *val++ );
5638 Real* vec, int* idx,
5639 Real* rhs, int* ridx, int rn,
5640 Real* forest, int* forestNum, int* forestIdx )
5652 int* it = forestIdx;
5656 for ( i = j = 0; i < rn; ++i )
5659 assert( k >= 0 && k < thedim );
5671 *forestNum = rn = j;
5681 for ( i = j = 0; i < rn; ++i )
5684 assert( k >= 0 && k < thedim );
5705 Real* vec, int* idx,
5706 Real* rhs, int* ridx, int rn,
5708 Real* rhs2, int* ridx2, int rn2,
5709 Real* forest, int* forestNum, int* forestIdx )
5715 vSolveLright2( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2 );
5725 int* it = forestIdx;
5729 for ( i = j = 0; i < rn; ++i )
5732 assert( k >= 0 && k < thedim );
5744 *forestNum = rn = j;
5754 for ( i = j = 0; i < rn; ++i )
5757 assert( k >= 0 && k < thedim );
5771 ridx2[0] = thedim - 1;
5784 for ( i = j = 0; i < rn2; ++i )
5787 assert( k >= 0 && k < thedim );
5828 Real* rhs, int* ridx, int& rn,
5830 Real* rhs2, int* ridx2, int& rn2,
5831 Real* forest, int* forestNum, int* forestIdx )
5834 vSolveLright2( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2 );
5843 int* it = forestIdx;
5845 for ( i = j = 0; i < rn; ++i )
5848 assert( k >= 0 && k < thedim );
5860 *forestNum = rn = j;
5864 for ( i = j = 0; i < rn; ++i )
5867 assert( k >= 0 && k < thedim );
5879 for ( i = j = 0; i < rn2; ++i )
5882 assert( k >= 0 && k < thedim );
5895 rn2 = vSolveUright( vec2, idx2, rhs2, ridx2, rn2, eps2 );
5906 Real* vec, int* idx,
5907 Real* rhs, int* ridx, int rn,
5909 Real* rhs2, int* ridx2, int rn2,
5911 Real* rhs3, int* ridx3, int rn3,
5912 Real* forest, int* forestNum, int* forestIdx )
5915 vSolveLright3( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2, rhs3, ridx3, &rn3, eps3 );
5916 assert( rn >= 0 && rn <= thedim );
5917 assert( rn2 >= 0 && rn2 <= thedim );
5918 assert( rn3 >= 0 && rn3 <= thedim );
5928 int* it = forestIdx;
5932 for ( i = j = 0; i < rn; ++i )
5935 assert( k >= 0 && k < thedim );
5947 *forestNum = rn = j;
5957 for ( i = j = 0; i < rn; ++i )
5960 assert( k >= 0 && k < thedim );
5974 ridx2[0] = thedim - 1;
5984 for ( i = j = 0; i < rn2; ++i )
5987 assert( k >= 0 && k < thedim );
6006 if ( rn3 > thedim*verySparseFactor4right )
6008 ridx3[0] = thedim - 1;
6018 for ( i = j = 0; i < rn3; ++i )
6021 assert( k >= 0 && k < thedim );
6056 Real* rhs, int* ridx, int& rn,
6058 Real* rhs2, int* ridx2, int& rn2,
6060 Real* rhs3, int* ridx3, int& rn3,
6061 Real* forest, int* forestNum, int* forestIdx )
6063 vSolveLright3( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2, rhs3, ridx3, &rn3, eps3 );
6064 assert( rn >= 0 && rn <= thedim );
6065 assert( rn2 >= 0 && rn2 <= thedim );
6066 assert( rn3 >= 0 && rn3 <= thedim );
6075 int* it = forestIdx;
6077 for ( i = j = 0; i < rn; ++i )
6080 assert( k >= 0 && k < thedim );
6092 *forestNum = rn = j;
6096 for ( i = j = 0; i < rn; ++i )
6099 assert( k >= 0 && k < thedim );
6111 for ( i = j = 0; i < rn2; ++i )
6114 assert( k >= 0 && k < thedim );
6125 for ( i = j = 0; i < rn3; ++i )
6128 assert( k >= 0 && k < thedim );
6140 rn2 = vSolveUright( vec2, idx2, rhs2, ridx2, rn2, eps2 );
6141 rn3 = vSolveUright( vec3, idx3, rhs3, ridx3, rn3, eps3 );
6153 Real* rhs2, int* ridx2, int rn2 )
6156 assert( rn2 >= 0 && rn2 <= thedim );
6160 *ridx2 = thedim - 1;
6172 for ( i = j = 0; i < rn2; ++i )
6175 assert( k >= 0 && k < thedim );
6205 Real* vec, int* idx,
6206 Real* rhs, int* ridx, int rn )
6212 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6216 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6233 Real* vec, int* idx,
6234 Real* rhs, int* ridx, int rn,
6236 Real* rhs2, int* ridx2, int rn2 )
6242 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6248 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6262 Real* vec, int* idx,
6263 Real* rhs, int* ridx, int& rn,
6264 Real* vec2, int* idx2,
6265 Real* rhs2, int* ridx2, int& rn2 )
6270 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6272 rn2 = solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6276 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6278 rn2 = solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6289 Real* vec, int* idx,
6290 Real* rhs, int* ridx, int rn,
6292 Real* rhs2, int* ridx2, int rn2,
6294 Real* rhs3, int* ridx3, int rn3 )
6300 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6308 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6325 Real* vec, int* idx,
6326 Real* rhs, int* ridx, int& rn,
6327 Real* vec2, int* idx2,
6328 Real* rhs2, int* ridx2, int& rn2,
6329 Real* vec3, int* idx3,
6330 Real* rhs3, int* ridx3, int& rn3 )
6335 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6337 rn2 = solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6339 rn3 = solveUleft( eps, vec3, idx3, rhs3, ridx3, rn3 );
6343 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6345 rn2 = solveUleft( eps, vec2, idx2, rhs2, ridx2, rn2 );
6347 rn3 = solveUleft( eps, vec3, idx3, rhs3, ridx3, rn3 );
6359 Real* rhs2, int* ridx2, int rn2 )
|