|
SoPlex Doxygen Documentation
|
Go to the documentation of this file.
60 for ( i = 1; i < *size; ++i )
61 for ( j = 0; j < i; ++j )
62 assert( heap[i] != heap[j] );
74 e = heap[s = --( *size )];
77 for ( j = 0, i = 1; i < s; i = 2 * j + 1 )
110 if ( i < *size && e < heap[i] )
131 if ( elem < heap[i] )
144 for ( i = 1; i < *size; ++i )
145 for ( j = 0; j < i; ++j )
146 assert( heap[i] != heap[j] );
158 e = heap[s = --( *size )];
161 for ( j = 0, i = 1; i < s; i = 2 * j + 1 )
194 if ( i < *size && e > heap[i] )
236 if ( pivot_col != 0 )
239 if ( pivot_colNZ != 0 )
242 if ( pivot_row != 0 )
245 if ( pivot_rowNZ != 0 )
257 METHOD( "CLUFactor::initPerm()" );
259 for ( int i = 0; i < thedim; ++i )
270 METHOD( "CLUFactor::setPivot()" );
290 METHOD( "CLUFactor::packRows()" );
303 for ( ring = list-> next; ring != list; ring = ring-> next )
307 if ( l_rbeg[l_row] != n )
313 assert( l_rlen[l_row] <= l_rmax[l_row] );
315 l_rmax[l_row] = l_rlen[l_row];
316 j = i + l_rlen[l_row];
318 for ( ; i < j; ++i, ++n )
321 l_ridx[n] = l_ridx[i];
322 l_rval[n] = l_rval[i];
327 while ( ring != list );
329 goto terminatePackRows;
334 l_rmax[l_row] = l_rlen[l_row];
349 METHOD( "CLUFactor::forestPackColumns()" );
362 for ( ring = list-> next; ring != list; ring = ring-> next )
366 if ( cbeg[colno] != n )
373 cmax[colno] = clen[colno];
384 while ( ring != list );
386 goto terminatePackColumns;
391 cmax[colno] = clen[colno];
394 terminatePackColumns :
405 METHOD( "CLUFactor::remaxRow()" );
406 assert( u. row. max[p_row] < len );
410 int delta = len - u. row. max[p_row];
415 delta = len - u. row. max[p_row];
425 && "ERROR: could not allocate memory for row file" );
449 && "ERROR: could not allocate memory for row file" );
466 for ( ; i < k; ++i, ++j )
484 METHOD( "CLUFactor::packColumns()" );
496 for ( ring = list-> next; ring != list; ring = ring-> next )
500 if ( l_cbeg[l_col] != n )
507 l_cmax[l_col] = l_clen[l_col];
508 j = i + l_clen[l_col];
511 l_cidx[n++] = l_cidx[i];
515 while ( ring != list );
517 goto terminatePackColumns;
522 l_cmax[l_col] = l_clen[l_col];
525 terminatePackColumns :
536 METHOD( "CLUFactor::remaxCol()" );
537 assert( u. col. max[p_col] < len );
541 int delta = len - u. col. max[p_col];
546 delta = len - u. col. max[p_col];
556 && "ERROR: could not allocate memory for column file" );
580 && "ERROR: could not allocate memory for column file" );
606 METHOD( "CLUFactor::forestReMaxCol()" );
607 assert( u. col. max[p_col] < len );
611 int delta = len - u. col. max[p_col];
616 delta = len - u. col. max[p_col];
624 && "ERROR: could not allocate memory for column file" );
647 && "ERROR: could not allocate memory for column file" );
706 METHOD( "CLUFactor::forestUpdate()" );
707 int i, j, k, h, m, n;
741 for ( i += j - 1; i >= j; --i )
745 h = --( rlen[m] ) + k;
747 while ( ridx[k] != p_col )
772 if ( num > cmax[p_col] )
783 for ( j = 0; j < num; ++j )
791 if ( fabs( x ) > l_maxabs )
792 l_maxabs = fabs( x );
795 assert( k - cbeg[p_col] < cmax[p_col] );
802 if ( rmax[i] <= rlen[i] )
809 h = rbeg[i] + ( rlen[i] )++;
821 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
838 for ( i = 0; i < dim; ++i )
845 if ( fabs( x ) > l_maxabs )
846 l_maxabs = fabs( x );
851 clen[p_col] = k - cbeg[p_col];
860 assert( k - cbeg[p_col] < cmax[p_col] );
867 if ( rmax[i] <= rlen[i] )
874 h = rbeg[i] + ( rlen[i] )++;
886 nzCnt += ( clen[p_col] = k - cbeg[p_col] );
888 if ( cbeg[p_col] + cmax[p_col] == u. col. used )
891 cmax[p_col] = clen[p_col];
904 for ( i = c; i < r; ++i )
905 rorig[i] = rorig[i + 1];
909 for ( i = c; i <= r; ++i )
914 for ( i = c; i < r; ++i )
915 corig[i] = corig[i + 1];
919 for ( i = c; i <= r; ++i )
948 for ( i += j - 1; i >= j; --i )
953 m = --( clen[k] ) + cbeg[k];
955 for ( h = m; cidx[h] != rowno; --h )
972 assert(( num == 0 ) || ( nonz != 0 ) );
980 for ( i = 0; i < num; ++i )
981 assert( p_work[corig[nonz[i]]] != 0.0 );
991 assert( p_work[k] != 0.0 );
995 x = p_work[k] * diag[n];
1005 if ( fabs( x ) > l_maxabs )
1006 l_maxabs = fabs( x );
1012 for ( ; j < m; ++j )
1015 Real y = p_work[jj];
1022 p_work[jj] = y + (( y == 0 ) ? MARKER : 0.0 );
1043 diag[rowno] = 1 / x;
1050 if ( rmax[rowno] < num )
1064 for ( i = 0; i < num; ++i )
1077 if ( fabs( x ) > l_maxabs )
1078 l_maxabs = fabs( x );
1088 if ( clen[j] >= cmax[j] )
1095 cval[cbeg[j] + clen[j]] = x;
1097 cidx[cbeg[j] + clen[j]++] = rowno;
1101 rlen[rowno] = n - rbeg[rowno];
1107 for ( i += j - 1; i >= j; --i )
1110 p_work[k] = rval[i];
1111 m = --( clen[k] ) + cbeg[k];
1113 for ( h = m; cidx[h] != rowno; --h )
1130 for ( i = c; i < r; ++i )
1137 x = p_work[k] * diag[n];
1143 if ( fabs( x ) > l_maxabs )
1144 l_maxabs = fabs( x );
1150 for ( ; j < m; ++j )
1151 p_work[ridx[j]] -= x * rval[j];
1174 diag[rowno] = 1 / x;
1184 for ( i = r + 1; i < dim; ++i )
1185 if ( p_work[corig[i]] != 0.0 )
1188 if ( rmax[rowno] < n )
1202 for ( i = r + 1; i < dim; ++i )
1209 if ( fabs( x ) > l_maxabs )
1210 l_maxabs = fabs( x );
1220 if ( clen[j] >= cmax[j] )
1227 cval[cbeg[j] + clen[j]] = x;
1229 cidx[cbeg[j] + clen[j]++] = rowno;
1233 rlen[rowno] = n - rbeg[rowno];
1244 i = rbeg[rowno] + --( rlen[rowno] );
1245 diag[rowno] = 1 / rval[i];
1247 for ( j = i = --( clen[p_col] ) + cbeg[p_col]; cidx[i] != rowno; --i )
1269 METHOD( "CLUFactor::update()" );
1275 assert( p_work[p_col] != 0 );
1276 rezi = 1 / p_work[p_col];
1284 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1287 lval[ll] = rezi * p_work[j];
1294 lval[ll] = 1 - rezi;
1297 for ( --i; i >= 0; --i )
1301 lval[ll] = x = rezi * p_work[j];
1305 if ( fabs( x ) > maxabs )
1318 METHOD( "CLUFactor::updateNoClear()" );
1324 assert( p_work[p_col] != 0 );
1325 rezi = 1 / p_work[p_col];
1331 for ( i = num - 1; ( j = p_idx[i] ) != p_col; --i )
1334 lval[ll] = rezi * p_work[j];
1340 lval[ll] = 1 - rezi;
1343 for ( --i; i >= 0; --i )
1347 lval[ll] = x = rezi * p_work[j];
1350 if ( fabs( x ) > maxabs )
1385 METHOD( "CLUFactor::initFactorMatrix()" );
1390 Dring *rring, *lastrring;
1391 Dring *cring, *lastcring;
1401 for ( int i = 0; i < thedim; i++ )
1406 for ( int i = 0; i < thedim; i++ )
1417 for ( int j = 0; j < k; ++j )
1450 lastrring-> next = rring;
1458 lastcring-> next = cring;
1462 for ( int i = 0; i < thedim; i++ )
1468 rring-> prev = lastrring;
1469 lastrring-> next = rring;
1474 cring-> prev = lastcring;
1475 lastcring-> next = cring;
1499 for ( int i = 0; i < thedim; i++ )
1509 for ( int j = 0; j < psv-> size() && nnonzeros <= 1; j++ )
1516 if ( nnonzeros == 0 )
1524 if ( nnonzeros == 1 )
1530 for ( j = 0; isZero( psv-> value( j ), eps ); j++ )
1533 assert( j < psv->size() );
1543 x = psv-> value( j );
1566 assert( nnonzeros >= 2 );
1569 for ( int j = 0; j < psv-> size(); j++ )
1571 x = psv-> value( j );
1576 k = psv-> index( j );
1595 assert( nnonzeros >= 2 );
1614 METHOD( "CLUFactor::colSingletons()" );
1617 int p_col, p_row, newrow;
1633 assert( p_row >= 0 );
1637 for ( j = 0; j < len; ++j )
1646 for ( k = n; u. col. idx[k] != p_row; ++k )
1664 if ( rperm[newrow] >= 0 )
1674 for ( k = n; u. row. idx[k] != p_col; --k )
1708 METHOD( "CLUFactor::rowSingletons()" );
1711 int p_row, p_col, len, rs, lk;
1720 for ( i = 0; i < thedim; ++i )
1722 if ( rperm[i] < 0 && u. row. len[i] == 1 )
1738 setPivot( rs, p_col, p_row, pval );
1748 i = ( u. col. len[p_col] -= i );
1750 for ( ; i < len; ++i )
1761 for ( j = k; u. row. idx[j] != p_col; --j )
1804 METHOD( "CLUFactor::initFactorRings()" );
1822 for ( i = 0; i < thedim; ++i )
1858 METHOD( "CLUFactor::freeFactorRings()" );
1880 METHOD( "CLUFactor::eliminateRowSingletons()" );
1909 for ( ; ( r = idx[i] ) != prow; ++i )
1916 for ( j = k; u. row. idx[j] != pcol; --j )
1947 assert( i < len && "ERROR: pivot column does not contain pivot row" );
1949 for ( ++i; i < len; ++i )
1957 for ( j = k; u. row. idx[j] != pcol; --j )
2002 METHOD( "CLUFactor::eliminateColSingletons()" );
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 )
2080 METHOD( "CLUFactor::selectPivots()" );
2101 if ( candidates > 4 )
2116 len = u. row. len[rw] + beg - 1;
2124 l_maxabs = fabs( u. row. val[len] );
2126 for ( i = len - 1; i >= beg; --i )
2127 if ( l_maxabs < fabs( u. row. val[i] ) )
2128 l_maxabs = fabs( u. row. val[i] );
2133 l_maxabs *= threshold;
2139 for ( i = len; i >= beg; --i )
2145 if ( j < mkwtz && fabs( x ) > l_maxabs )
2161 len = u. col. len[cl] + beg - 1;
2169 for ( i = len; i >= beg; --i )
2200 l_maxabs = fabs( u. row. val[m] );
2202 for ( kk = m + u. row. len[k] - 1; kk >= m; --kk )
2204 if ( l_maxabs < fabs( u. row. val[kk] ) )
2205 l_maxabs = fabs( u. row. val[kk] );
2215 for ( --kk; kk > m; --kk )
2217 if ( l_maxabs < fabs( u. row. val[kk] ) )
2218 l_maxabs = fabs( u. row. val[kk] );
2224 l_maxabs *= threshold;
2226 if ( fabs( x ) > l_maxabs )
2232 if ( j <= count + 1 )
2261 if ( pr-> idx == rw || pr-> mkwtz >= mkwtz )
2267 if ( pr-> idx != rw )
2275 if ( num >= candidates )
2306 METHOD( "CLUFactor::updateRow()" );
2309 int c, i, j, k, ll, m, n;
2318 for ( j = m; u. row. idx[j] != pcol; --j )
2339 for ( j = m - 1; j >= n; --j )
2368 for ( i = k; u. col. idx[i] != r; --i )
2381 if ( ll + fill > u. row. max[r] )
2438 METHOD( "CLUFactor::eliminatePivot()" );
2439 int i, j, k, m = -1;
2444 int plen = --( u. row. len[prow] );
2445 int pend = pbeg + plen;
2474 for ( i = pbeg; i < pend; ++i )
2482 for ( k = m; u. col. idx[k] != prow; ++k )
2500 updateRow( m, lv++, prow, pcol, pval, eps );
2507 for ( ++i; i < m; ++i )
2519 for ( i = u. row. start[prow], pend = i + plen; i < pend; ++i )
2534 const Real threshold )
2536 METHOD( "CLUFactor::eliminateNucleus()" );
2555 for ( i = 0; i < thedim; ++i )
2574 "ERROR: no pivot element selected" );
2577 pivot = pivot-> next )
2597 "ERROR: one row must be left" );
2599 "ERROR: one col must be left" );
2612 METHOD( "CLUFactor::setupColVals()" );
2622 for ( i = 0; i < thedim; i++ )
2627 for ( i = 0; i < thedim; i++ )
2638 assert(( *idx >= 0 ) && ( *idx < thedim ) );
2642 assert(( k >= 0 ) && ( k < u. col. size ) );
2651 if ( fabs( *val ) > maxabs )
2668 METHOD( "CLUFactor::setupRowVals()" );
2729 for ( i = thedim; i--; *l_rbeg++ = 0 )
2731 *rorig++ = *rrorig++;
2732 *rperm++ = *rrperm++;
2737 l_rbeg = l. rbeg + 1;
2739 for ( i = mem; i--; )
2744 for ( m = 0, i = thedim; i--; l_rbeg++ )
2753 l_rbeg = l. rbeg + 1;
2755 for ( i = j = 0; i < vecs; ++i )
2760 for ( ; j < beg[i + 1]; j++ )
2762 k = l_rbeg[*idx++]++;
2771 assert( l. rbeg[0] == 0 );
2782 METHOD( "CLUFactor::factor()" );
2838 METHOD( "CLUFactor::dump()" );
2849 for ( i = 0; i < thedim; ++i )
2853 << "] = " << diag[i] << std::endl;
2855 for ( j = 0; j < u. row. len[i]; ++j )
2856 spxout << "DCLUFA02 u[" << i << "]: ["
2863 for ( i = 0; i < thedim; ++i )
2868 spxout << "DCLUFA03 l[" << i << "]" << std::endl;
2872 << "]: [" << l. idx[k]
2873 << "] = " << l. val[k] << std::endl;
2890 METHOD( "CLUFactor::minRowMem()" );
2906 METHOD( "CLUFactor::minColMem()" );
2917 METHOD( "CLUFactor::forestMinColMem()" );
2929 METHOD( "CLUFactor::minLMem()" );
2931 if ( size > l. size )
2942 METHOD( "CLUFactor::makeLvec()" );
2950 int* p_lrow = l. row;
2955 assert( p_len > 0 && "ERROR: no empty columns allowed in L vectors" );
2971 #ifdef ENABLE_CONSISTENCY_CHECKS
2972 METHOD( "CLUFactor::isConsistent()" );
2989 assert( ring-> idx >= 0 );
2991 assert( ring-> prev-> next == ring );
3014 assert( ring-> idx >= 0 );
3016 assert( ring-> prev-> next == ring );
3031 for ( i = 0; i < thedim; ++i )
3040 for ( i = 0; i < thedim; ++i )
3067 for ( i = 0; i < thedim; ++i )
3089 for ( i = 0; i < thedim - temp. stage; ++i )
3103 #if 0 // Stimmt leider nicht
3108 for ( i = 0; i < thedim; i++ )
3110 assert( l. rbeg[i] >= 0 );
3114 assert( l. rorig[i] >= 0 );
3115 assert( l. rorig[i] < thedim );
3116 assert( l. rperm[i] >= 0 );
3117 assert( l. rperm[i] < thedim );
3119 assert( l. ridx[i] >= 0 );
3120 assert( l. ridx[i] < thedim );
3125 #endif // CONSISTENCY_CHECKS
3132 METHOD( "CLUFactor::solveUright()" );
3134 for ( int i = thedim - 1; i >= 0; i-- )
3138 Real x = wrk[c] = diag[r] * vec[r];
3153 METHOD( "CLUFactor::solveUrightEps()" );
3156 int *cidx, *clen, *cbeg;
3173 for ( i = thedim - 1; i >= 0; --i )
3176 x = diag[r] * rhs[r];
3183 val = &cval[cbeg[c]];
3184 idx = &cidx[cbeg[c]];
3188 rhs[*idx++] -= x * ( *val++ );
3198 METHOD( "CLUFactor::solveUright2()" );
3201 int *cidx, *clen, *cbeg;
3216 for ( i = thedim - 1; i >= 0; --i )
3220 p_work1[c] = x1 = diag[r] * vec1[r];
3221 p_work2[c] = x2 = diag[r] * vec2[r];
3222 vec1[r] = vec2[r] = 0;
3224 if ( x1 != 0.0 && x2 != 0.0 )
3226 val = &cval[cbeg[c]];
3227 idx = &cidx[cbeg[c]];
3232 vec1[*idx] -= x1 * ( *val );
3233 vec2[*idx++] -= x2 * ( *val++ );
3239 val = &cval[cbeg[c]];
3240 idx = &cidx[cbeg[c]];
3244 vec1[*idx++] -= x1 * ( *val++ );
3249 val = &cval[cbeg[c]];
3250 idx = &cidx[cbeg[c]];
3254 vec2[*idx++] -= x2 * ( *val++ );
3261 int* nonz, Real eps )
3263 METHOD( "CLUFactor::solveUright2eps()" );
3266 int *cidx, *clen, *cbeg;
3267 bool notzero1, notzero2;
3284 for ( i = thedim - 1; i >= 0; --i )
3288 p_work1[c] = x1 = diag[r] * vec1[r];
3289 p_work2[c] = x2 = diag[r] * vec2[r];
3290 vec1[r] = vec2[r] = 0;
3294 if ( notzero1 && notzero2 )
3298 val = &cval[cbeg[c]];
3299 idx = &cidx[cbeg[c]];
3304 vec1[*idx] -= x1 * ( *val );
3305 vec2[*idx++] -= x2 * ( *val++ );
3314 val = &cval[cbeg[c]];
3315 idx = &cidx[cbeg[c]];
3319 vec1[*idx++] -= x1 * ( *val++ );
3325 val = &cval[cbeg[c]];
3326 idx = &cidx[cbeg[c]];
3330 vec2[*idx++] -= x2 * ( *val++ );
3344 METHOD( "CLUFactor::solveLright()" );
3349 int *lrow, *lidx, *idx;
3359 for ( i = 0; i < end; ++i )
3361 if (( x = vec[lrow[i]] ) != 0.0 )
3363 MSG_DEBUG( spxout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl; )
3369 for ( j = lbeg[i + 1]; j > k; --j )
3371 MSG_DEBUG( spxout << " -> y" << *idx << " -= " << x << " * " << *val << " = " << x * ( *val ) << " -> " << vec[*idx] - x * ( *val ) << std::endl; )
3372 vec[*idx++] -= x * ( *val++ );
3383 for ( ; i < end; ++i )
3390 for ( j = lbeg[i + 1]; j > k; --j )
3391 x += vec[*idx++] * ( *val++ );
3395 MSG_DEBUG( spxout << "y" << lrow[i] << "=" << vec[lrow[i]] << std::endl; )
3404 METHOD( "CLUFactor::solveLright2()" );
3410 int *lrow, *lidx, *idx;
3420 for ( i = 0; i < end; ++i )
3425 if ( x1 != 0 && x2 != 0.0 )
3431 for ( j = lbeg[i + 1]; j > k; --j )
3433 vec1[*idx] -= x1 * ( *val );
3434 vec2[*idx++] -= x2 * ( *val++ );
3444 for ( j = lbeg[i + 1]; j > k; --j )
3445 vec1[*idx++] -= x1 * ( *val++ );
3454 for ( j = lbeg[i + 1]; j > k; --j )
3455 vec2[*idx++] -= x2 * ( *val++ );
3463 for ( ; i < end; ++i )
3471 for ( j = lbeg[i + 1]; j > k; --j )
3473 x1 += vec1[*idx] * ( *val );
3474 x2 += vec2[*idx++] * ( *val++ );
3477 vec1[lrow[i]] -= x1;
3479 vec2[lrow[i]] -= x2;
3486 METHOD( "CLUFactor::solveUpdateRight()" );
3491 int *lrow, *lidx, *idx;
3505 if (( x = vec[lrow[i]] ) != 0.0 )
3511 for ( j = lbeg[i + 1]; j > k; --j )
3512 vec[*idx++] -= x * ( *val++ );
3519 METHOD( "CLUFactor::solveUpdateRight2()" );
3544 if ( x1 != 0.0 && x2 != 0.0 )
3550 for ( j = lbeg[i + 1]; j > k; --j )
3552 vec1[*idx] -= x1 * ( *val );
3553 vec2[*idx++] -= x2 * ( *val++ );
3563 for ( j = lbeg[i + 1]; j > k; --j )
3564 vec1[*idx++] -= x1 * ( *val++ );
3573 for ( j = lbeg[i + 1]; j > k; --j )
3574 vec2[*idx++] -= x2 * ( *val++ );
3580 Real* rhs, Real* forest, int* forestNum, int* forestIdx )
3582 METHOD( "CLUFactor::solveRight4update()" );
3589 for ( int i = 0; i < thedim; i++ )
3593 n += rhs[i] != 0.0 ? 1 : 0;
3611 METHOD( "CLUFactor::solveRight()" );
3629 METHOD( "CLUFactor::solveRight2update()" );
3636 for ( int i = 0; i < thedim; i++ )
3639 forest[i] = rhs1[i];
3640 n += rhs1[i] != 0.0 ? 1 : 0;
3662 METHOD( "CLUFactor::solveRight2()" );
3678 METHOD( "CLUFactor::solveUleft()" );
3682 int *ridx, *rlen, *rbeg, *idx;
3693 for ( i = 0; i < thedim; ++i )
3708 for ( int m = rlen[r]; m != 0; --m )
3709 vec[*idx++] -= x * ( *val++ );
3717 METHOD( "CLUFactor::solveUleft()" );
3720 for ( int i = 0; i < thedim; ++i )
3725 for ( int m = u. row. start[r]; m < end; m++ )
3727 std::cout << u. row. val[m] << " ";
3734 for ( int i = 0; i < thedim; ++i )
3747 #if defined(WITH_WARNINGS) || !defined(NDEBUG)
3764 for ( int m = u. row. start[r]; m < end; m++ )
3782 METHOD( "CLUFactor::solveUleft2()" );
3787 int *ridx, *rlen, *rbeg, *idx;
3798 for ( i = 0; i < thedim; ++i )
3805 if (( x1 != 0.0 ) && ( x2 != 0.0 ) )
3815 for ( int m = rlen[r]; m != 0; --m )
3817 vec1[*idx] -= x1 * ( *val );
3818 vec2[*idx] -= x2 * ( *val++ );
3831 for ( int m = rlen[r]; m != 0; --m )
3832 vec1[*idx++] -= x1 * ( *val++ );
3843 for ( int m = rlen[r]; m != 0; --m )
3844 vec2[*idx++] -= x2 * ( *val++ );
3855 METHOD( "CLUFactor::solveLleft2forest()" );
3862 int *lidx, *idx, *lrow;
3886 for ( j = lbeg[i + 1]; j > k; --j )
3888 vec1[*idx] -= x1 * ( *val );
3889 vec2[*idx++] -= x2 * ( *val++ );
3898 for ( j = lbeg[i + 1]; j > k; --j )
3899 vec1[*idx++] -= x1 * ( *val++ );
3909 for ( j = lbeg[i + 1]; j > k; --j )
3910 vec2[*idx++] -= x2 * ( *val++ );
3923 METHOD( "CLUFactor::solveLleft2()" );
3946 for ( ; i >= 0; --i )
3954 for ( j = lbeg[i + 1]; j > k; --j )
3956 x1 += vec1[*idx] * ( *val );
3957 x2 += vec2[*idx++] * ( *val++ );
3960 vec1[lrow[i]] -= x1;
3962 vec2[lrow[i]] -= x2;
3966 for ( i = thedim; i--; )
3971 x1not0 = ( x1 != 0 );
3972 x2not0 = ( x2 != 0 );
3974 if ( x1not0 && x2not0 )
3977 j = rbeg[r + 1] - k;
3984 vec1[*idx] -= x1 * *val;
3985 vec2[*idx++] -= x2 * *val++;
3992 j = rbeg[r + 1] - k;
3999 vec1[*idx++] -= x1 * *val++;
4006 j = rbeg[r + 1] - k;
4013 vec2[*idx++] -= x2 * *val++;
4023 METHOD( "CLUFactor::solveLleftForest()" );
4027 int *idx, *lidx, *lrow, *lbeg;
4038 if (( x = vec[lrow[i]] ) != 0.0 )
4044 for ( j = lbeg[i + 1]; j > k; --j )
4045 vec[*idx++] -= x * ( *val++ );
4054 METHOD( "CLUFactor::solveLleft()" );
4071 for ( int j = lbeg[i + 1]; j > k; --j )
4072 x += vec[*idx++] * ( *val++ );
4078 for ( int i = thedim - 1; i >= 0; --i )
4085 for ( int k = l. rbeg[r]; k < l. rbeg[r + 1]; k++ )
4089 assert( l. rperm[j] < i );
4091 vec[j] -= x * l. rval[k];
4096 #endif // WITH_L_ROWS
4101 METHOD( "CLUFactor::solveLleftEps()" );
4128 for ( j = lbeg[i + 1]; j > k; --j )
4129 x += vec[*idx++] * ( *val++ );
4135 for ( i = thedim; i--; )
4145 j = rbeg[r + 1] - k;
4152 vec[*idx++] -= x * *val++;
4166 METHOD( "CLUFactor::solveUpdateLeft()" );
4170 int *lrow, *lidx, *idx;
4189 for ( j = lbeg[i + 1]; j > k; --j )
4190 x += vec[*idx++] * ( *val++ );
4198 METHOD( "CLUFactor::solveUpdateLeft2()" );
4202 int *lrow, *lidx, *idx;
4222 for ( j = lbeg[i + 1]; j > k; --j )
4224 x1 += vec1[*idx] * ( *val );
4225 x2 += vec2[*idx++] * ( *val++ );
4228 vec1[lrow[i]] -= x1;
4230 vec2[lrow[i]] -= x2;
4236 METHOD( "CLUFactor::solveUpdateLeft()" );
4240 int *lrow, *lidx, *idx;
4255 assert( k >= 0 && k < l. size );
4260 for ( j = lbeg[i + 1]; j > k; --j )
4262 assert( *idx >= 0 && *idx < thedim );
4263 x += vec[*idx++] * ( *val++ );
4283 vec[k] = ( y != 0 ) ? y : MARKER;
4292 METHOD( "CLUFactor::solveLeft()" );
4311 METHOD( "CLUFactor::solveLeftEps()" );
4335 METHOD( "CLUFactor::solveLeft2()" );
4354 Real* vec, int* vecidx,
4355 Real* rhs, int* rhsidx, int rhsn )
4357 METHOD( "CLUFactor::solveUleft()" );
4359 int i, j, k, n, r, c;
4360 int *rorig, *corig, *cperm;
4361 int *ridx, *rlen, *rbeg, *idx;
4371 for ( i = 0; i < rhsn; )
4387 assert( i >= 0 && i < thedim );
4389 assert( c >= 0 && c < thedim );
4396 assert( r >= 0 && r < thedim );
4401 assert( k >= 0 && k < u. row. size );
4405 for ( int m = rlen[r]; m; --m )
4408 assert( j >= 0 && j < thedim );
4413 y = -x * ( *val++ );
4423 y -= x * ( *val++ );
4424 rhs[j] = ( y != 0 ) ? y : MARKER;
4435 Real* rhs, int* rhsidx, int rhsn )
4437 METHOD( "CLUFactor::solveUleftNoNZ()" );
4440 int *rorig, *corig, *cperm;
4441 int *ridx, *rlen, *rbeg, *idx;
4451 for ( i = 0; i < rhsn; )
4465 assert( i >= 0 && i < thedim );
4467 assert( c >= 0 && c < thedim );
4474 assert( r >= 0 && r < thedim );
4478 assert( k >= 0 && k < u. row. size );
4482 for ( int m = rlen[r]; m; --m )
4485 assert( j >= 0 && j < thedim );
4490 y = -x * ( *val++ );
4500 y -= x * ( *val++ );
4501 rhs[j] = ( y != 0 ) ? y : MARKER;
4511 METHOD( "CLUFactor::solveLleftForest()" );
4515 int *idx, *lidx, *lrow, *lbeg;
4525 assert( i >= 0 && i < l. size );
4527 if (( x = vec[lrow[i]] ) != 0.0 )
4530 assert( k >= 0 && k < l. size );
4534 for ( j = lbeg[i + 1]; j > k; --j )
4537 assert( m >= 0 && m < thedim );
4542 y = -x * ( *val++ );
4552 y -= x * ( *val++ );
4553 vec[m] = ( y != 0 ) ? y : MARKER;
4565 METHOD( "CLUFactor::solveLleftForestNoNZ()" );
4569 int *idx, *lidx, *lrow, *lbeg;
4579 if (( x = vec[lrow[i]] ) != 0.0 )
4581 assert( i >= 0 && i < l. size );
4583 assert( k >= 0 && k < l. size );
4587 for ( j = lbeg[i + 1]; j > k; --j )
4589 assert( *idx >= 0 && *idx < thedim );
4590 vec[*idx++] -= x * ( *val++ );
4599 METHOD( "CLUFactor::solveLleft()" );
4618 #pragma warn "Not yet implemented, define WITH_L_ROWS"
4624 for ( ; i >= 0; --i )
4631 for ( j = lbeg[i + 1]; j > k; --j )
4632 x += vec[*idx++] * ( *val++ );
4640 for ( i = 0; i < rn; )
4656 j = rbeg[r + 1] - k;
4662 assert( l. rperm[*idx] < i );
4679 vec[m] = ( y != 0 ) ? y : MARKER;
4687 for ( i = 0; i < n; ++i )
4698 METHOD( "CLUFactor::solveLleftNoNZ()" );
4719 assert( i < thedim );
4721 for ( ; i >= 0; --i )
4724 assert( k >= 0 && k < l. size );
4729 for ( j = lbeg[i + 1]; j > k; --j )
4731 assert( *idx >= 0 && *idx < thedim );
4732 x += vec[*idx++] * ( *val++ );
4739 for ( i = thedim; i--; )
4747 j = rbeg[r + 1] - k;
4753 assert( l. rperm[*idx] < i );
4754 vec[*idx++] -= x * *val++;
4764 METHOD( "CLUFactor::vSolveLright()" );
4769 int *lrow, *lidx, *idx;
4779 for ( i = 0; i < end; ++i )
4789 for ( j = lbeg[i + 1]; j > k; --j )
4791 assert( *idx >= 0 && *idx < thedim );
4792 ridx[rn] = n = *idx++;
4793 rn += ( vec[n] == 0 ) ? 1 : 0;
4794 vec[n] -= x * ( *val++ );
4795 vec[n] += ( vec[n] == 0 ) ? MARKER : 0;
4804 for ( ; i < end; ++i )
4811 for ( j = lbeg[i + 1]; j > k; --j )
4813 assert( *idx >= 0 && *idx < thedim );
4814 x += vec[*idx++] * ( *val++ );
4817 ridx[rn] = j = lrow[i];
4819 rn += ( vec[j] == 0 ) ? 1 : 0;
4821 vec[j] += ( vec[j] == 0 ) ? MARKER : 0;
4829 Real* vec, int* ridx, int* rnptr, Real eps,
4830 Real* vec2, int* ridx2, int* rn2ptr, Real eps2 )
4832 METHOD( "CLUFactor::vSolveLright2()" );
4838 int *lrow, *lidx, *idx;
4851 for ( i = 0; i < end; ++i )
4865 for ( j = lbeg[i + 1]; j > k; --j )
4867 assert( *idx >= 0 && *idx < thedim );
4868 ridx[rn] = ridx2[rn2] = n = *idx++;
4871 rn += ( y == 0 ) ? 1 : 0;
4872 rn2 += ( y2 == 0 ) ? 1 : 0;
4874 y2 -= x2 * ( *val++ );
4875 vec[n] = y + ( y == 0 ? MARKER : 0 );
4876 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
4881 for ( j = lbeg[i + 1]; j > k; --j )
4883 assert( *idx >= 0 && *idx < thedim );
4884 ridx[rn] = n = *idx++;
4886 rn += ( y == 0 ) ? 1 : 0;
4887 y -= x * ( *val++ );
4888 vec[n] = y + ( y == 0 ? MARKER : 0 );
4899 for ( j = lbeg[i + 1]; j > k; --j )
4901 assert( *idx >= 0 && *idx < thedim );
4902 ridx2[rn2] = n = *idx++;
4904 rn2 += ( y2 == 0 ) ? 1 : 0;
4905 y2 -= x2 * ( *val++ );
4906 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
4915 for ( ; i < end; ++i )
4922 for ( j = lbeg[i + 1]; j > k; --j )
4924 assert( *idx >= 0 && *idx < thedim );
4925 x += vec[*idx] * ( *val );
4926 x2 += vec2[*idx++] * ( *val++ );
4929 ridx[rn] = ridx2[rn2] = j = lrow[i];
4931 rn += ( vec[j] == 0 ) ? 1 : 0;
4932 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
4935 vec[j] += ( vec[j] == 0 ) ? MARKER : 0;
4936 vec2[j] += ( vec2[j] == 0 ) ? MARKER : 0;
4946 Real* vec, int* ridx, int* rnptr, Real eps,
4947 Real* vec2, int* ridx2, int* rn2ptr, Real eps2,
4948 Real* vec3, int* ridx3, int* rn3ptr, Real eps3 )
4950 METHOD( "CLUFactor::vSolveLright3()" );
4957 int *lrow, *lidx, *idx;
4971 for ( i = 0; i < end; ++i )
4989 for ( j = lbeg[i + 1]; j > k; --j )
4991 assert( *idx >= 0 && *idx < thedim );
4992 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++;
4996 rn += ( y == 0 ) ? 1 : 0;
4997 rn2 += ( y2 == 0 ) ? 1 : 0;
4998 rn3 += ( y3 == 0 ) ? 1 : 0;
5000 y2 -= x2 * ( *val );
5001 y3 -= x3 * ( *val++ );
5002 vec[n] = y + ( y == 0 ? MARKER : 0 );
5003 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
5004 vec3[n] = y3 + ( y3 == 0 ? MARKER : 0 );
5010 for ( j = lbeg[i + 1]; j > k; --j )
5012 assert( *idx >= 0 && *idx < thedim );
5013 ridx[rn] = ridx2[rn2] = n = *idx++;
5016 rn += ( y == 0 ) ? 1 : 0;
5017 rn2 += ( y2 == 0 ) ? 1 : 0;
5019 y2 -= x2 * ( *val++ );
5020 vec[n] = y + ( y == 0 ? MARKER : 0 );
5021 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
5029 for ( j = lbeg[i + 1]; j > k; --j )
5031 assert( *idx >= 0 && *idx < thedim );
5032 ridx[rn] = ridx3[rn3] = n = *idx++;
5035 rn += ( y == 0 ) ? 1 : 0;
5036 rn3 += ( y3 == 0 ) ? 1 : 0;
5038 y3 -= x3 * ( *val++ );
5039 vec[n] = y + ( y == 0 ? MARKER : 0 );
5040 vec3[n] = y3 + ( y3 == 0 ? MARKER : 0 );
5046 for ( j = lbeg[i + 1]; j > k; --j )
5048 assert( *idx >= 0 && *idx < thedim );
5049 ridx[rn] = n = *idx++;
5051 rn += ( y == 0 ) ? 1 : 0;
5052 y -= x * ( *val++ );
5053 vec[n] = y + ( y == 0 ? MARKER : 0 );
5067 for ( j = lbeg[i + 1]; j > k; --j )
5069 assert( *idx >= 0 && *idx < thedim );
5070 ridx2[rn2] = ridx3[rn3] = n = *idx++;
5073 rn2 += ( y2 == 0 ) ? 1 : 0;
5074 rn3 += ( y3 == 0 ) ? 1 : 0;
5075 y2 -= x2 * ( *val );
5076 y3 -= x3 * ( *val++ );
5077 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
5078 vec3[n] = y3 + ( y3 == 0 ? MARKER : 0 );
5084 for ( j = lbeg[i + 1]; j > k; --j )
5086 assert( *idx >= 0 && *idx < thedim );
5087 ridx2[rn2] = n = *idx++;
5089 rn2 += ( y2 == 0 ) ? 1 : 0;
5090 y2 -= x2 * ( *val++ );
5091 vec2[n] = y2 + ( y2 == 0 ? MARKER : 0 );
5103 for ( j = lbeg[i + 1]; j > k; --j )
5105 assert( *idx >= 0 && *idx < thedim );
5106 ridx3[rn3] = n = *idx++;
5108 rn3 += ( y3 == 0 ) ? 1 : 0;
5109 y3 -= x3 * ( *val++ );
5110 vec3[n] = y3 + ( y3 == 0 ? MARKER : 0 );
5119 for ( ; i < end; ++i )
5126 for ( j = lbeg[i + 1]; j > k; --j )
5128 assert( *idx >= 0 && *idx < thedim );
5129 x += vec[*idx] * ( *val );
5130 x2 += vec2[*idx] * ( *val );
5131 x3 += vec3[*idx++] * ( *val++ );
5134 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i];
5136 rn += ( vec[j] == 0 ) ? 1 : 0;
5137 rn2 += ( vec2[j] == 0 ) ? 1 : 0;
5138 rn3 += ( vec3[j] == 0 ) ? 1 : 0;
5142 vec[j] += ( vec[j] == 0 ) ? MARKER : 0;
5143 vec2[j] += ( vec2[j] == 0 ) ? MARKER : 0;
5144 vec3[j] += ( vec3[j] == 0 ) ? MARKER : 0;
5155 Real* rhs, int* ridx, int rn, Real eps )
5157 METHOD( "CLUFactor::vSolveUright()" );
5158 int i, j, k, r, c, n;
5161 int *cidx, *clen, *cbeg;
5184 assert( i >= 0 && i < thedim );
5186 assert( r >= 0 && r < thedim );
5188 x = diag[r] * rhs[r];
5194 assert( c >= 0 && c < thedim );
5197 val = &cval[cbeg[c]];
5198 idx = &cidx[cbeg[c]];
5203 assert( *idx >= 0 && *idx < thedim );
5205 assert( k >= 0 && k < thedim );
5210 y = -x * ( *val++ );
5220 y -= x * ( *val++ );
5221 y += ( y == 0 ) ? MARKER : 0;
5228 for ( i = *ridx; i >= 0; --i )
5231 assert( r >= 0 && r < thedim );
5232 x = diag[r] * rhs[r];
5238 assert( c >= 0 && c < thedim );
5241 val = &cval[cbeg[c]];
5242 idx = &cidx[cbeg[c]];
5247 assert( *idx >= 0 && *idx < thedim );
5248 rhs[*idx++] -= x * ( *val++ );
5263 Real* rhs, int* ridx, int rn, Real eps )
5265 METHOD( "CLUFactor::vSolveUrightNoNZ()" );
5269 int *cidx, *clen, *cbeg;
5289 for ( i = *ridx; i >= 0; --i )
5291 assert( i >= 0 && i < thedim );
5293 assert( r >= 0 && r < thedim );
5294 x = diag[r] * rhs[r];
5301 val = &cval[cbeg[c]];
5302 idx = &cidx[cbeg[c]];
5307 assert( *idx >= 0 && *idx < thedim );
5308 rhs[*idx++] -= x * ( *val++ );
5320 assert( i >= 0 && i < thedim );
5324 assert( r >= 0 && r < thedim );
5326 x = diag[r] * rhs[r];
5334 val = &cval[cbeg[c]];
5335 idx = &cidx[cbeg[c]];
5341 assert( k >= 0 && k < thedim );
5346 y = -x * ( *val++ );
5356 y -= x * ( *val++ );
5357 y += ( y == 0 ) ? MARKER : 0;
5367 Real* vec, int* vidx, Real* rhs, int* ridx, int rn, Real eps,
5368 Real* vec2, Real* rhs2, int* ridx2, int rn2, Real eps2 )
5370 METHOD( "CLUFactor::vSolveUright2()" );
5371 int i, j, k, r, c, n;
5374 int *cidx, *clen, *cbeg;
5393 while ( rn + rn2 > 0 )
5403 if ( *ridx2 > *ridx )
5406 if ( *ridx2 < *ridx )
5414 assert( i >= 0 && i < thedim );
5417 assert( r >= 0 && r < thedim );
5419 x = diag[r] * rhs[r];
5420 x2 = diag[r] * rhs2[r];
5430 val = &cval[cbeg[c]];
5431 idx = &cidx[cbeg[c]];
5439 assert( k >= 0 && k < thedim );
5444 y2 = -x2 * ( *val );
5454 y2 -= x2 * ( *val );
5455 rhs2[k] = ( y2 != 0 ) ? y2 : MARKER;
5462 y = -x * ( *val++ );
5472 y -= x * ( *val++ );
5473 y += ( y == 0 ) ? MARKER : 0;
5483 assert( k >= 0 && k < thedim );
5488 y = -x * ( *val++ );
5498 y -= x * ( *val++ );
5499 y += ( y == 0 ) ? MARKER : 0;
5509 assert( c >= 0 && c < thedim );
5511 val = &cval[cbeg[c]];
5512 idx = &cidx[cbeg[c]];
5518 assert( k >= 0 && k < thedim );
5523 y2 = -x2 * ( *val++ );
5533 y2 -= x2 * ( *val++ );
5534 rhs2[k] = ( y2 != 0 ) ? y2 : MARKER;
5541 if ( *ridx > *ridx2 )
5546 for ( ; i >= 0; --i )
5548 assert( i < thedim );
5550 assert( r >= 0 && r < thedim );
5551 x = diag[r] * rhs[r];
5552 x2 = diag[r] * rhs2[r];
5559 assert( c >= 0 && c < thedim );
5561 val = &cval[cbeg[c]];
5562 idx = &cidx[cbeg[c]];
5572 assert( *idx >= 0 && *idx < thedim );
5573 rhs[*idx] -= x * ( *val );
5574 rhs2[*idx++] -= x2 * ( *val++ );
5581 assert( *idx >= 0 && *idx < thedim );
5582 rhs2[*idx++] -= x2 * ( *val++ );
5590 assert( c >= 0 && c < thedim );
5593 val = &cval[cbeg[c]];
5594 idx = &cidx[cbeg[c]];
5599 assert( *idx >= 0 && *idx < thedim );
5600 rhs[*idx++] -= x * ( *val++ );
5614 METHOD( "CLUFactor::vSolveUpdateRight()" );
5619 int *lrow, *lidx, *idx;
5632 assert( i >= 0 && i < thedim );
5638 assert( k >= 0 && k < l. size );
5642 for ( j = lbeg[i + 1]; j > k; --j )
5644 int m = ridx[n] = *idx++;
5645 assert( m >= 0 && m < thedim );
5647 n += ( y == 0 ) ? 1 : 0;
5648 y = y - x * ( *val++ );
5649 vec[m] = ( y != 0 ) ? y : MARKER;
5659 METHOD( "CLUFactor::vSolveUpdateRightNoNZ()" );
5664 int *lrow, *lidx, *idx;
5677 assert( i >= 0 && i < thedim );
5679 if (( x = vec[lrow[i]] ) != 0.0 )
5682 assert( k >= 0 && k < l. size );
5686 for ( j = lbeg[i + 1]; j > k; --j )
5688 assert( *idx >= 0 && *idx < thedim );
5689 vec[*idx++] -= x * ( *val++ );
5697 Real* vec, int* idx,
5698 Real* rhs, int* ridx, int rn,
5699 Real* forest, int* forestNum, int* forestIdx )
5701 METHOD( "CLUFactor::vSolveRight4update()" );
5712 int* it = forestIdx;
5716 for ( i = j = 0; i < rn; ++i )
5719 assert( k >= 0 && k < thedim );
5731 *forestNum = rn = j;
5741 for ( i = j = 0; i < rn; ++i )
5744 assert( k >= 0 && k < thedim );
5765 Real* vec, int* idx,
5766 Real* rhs, int* ridx, int rn,
5768 Real* rhs2, int* ridx2, int rn2,
5769 Real* forest, int* forestNum, int* forestIdx )
5771 METHOD( "CLUFactor::vSolveRight4update2()" );
5776 vSolveLright2( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2 );
5786 int* it = forestIdx;
5790 for ( i = j = 0; i < rn; ++i )
5793 assert( k >= 0 && k < thedim );
5805 *forestNum = rn = j;
5815 for ( i = j = 0; i < rn; ++i )
5818 assert( k >= 0 && k < thedim );
5832 ridx2[0] = thedim - 1;
5845 for ( i = j = 0; i < rn2; ++i )
5848 assert( k >= 0 && k < thedim );
5889 Real* vec, int* idx,
5890 Real* rhs, int* ridx, int rn,
5892 Real* rhs2, int* ridx2, int rn2,
5894 Real* rhs3, int* ridx3, int rn3,
5895 Real* forest, int* forestNum, int* forestIdx )
5897 METHOD( "CLUFactor::vSolveRight4update3()" );
5899 vSolveLright3( rhs, ridx, &rn, eps, rhs2, ridx2, &rn2, eps2, rhs3, ridx3, &rn3, eps3 );
5900 assert( rn >= 0 && rn <= thedim );
5901 assert( rn2 >= 0 && rn2 <= thedim );
5902 assert( rn3 >= 0 && rn3 <= thedim );
5912 int* it = forestIdx;
5916 for ( i = j = 0; i < rn; ++i )
5919 assert( k >= 0 && k < thedim );
5931 *forestNum = rn = j;
5941 for ( i = j = 0; i < rn; ++i )
5944 assert( k >= 0 && k < thedim );
5958 ridx2[0] = thedim - 1;
5968 for ( i = j = 0; i < rn2; ++i )
5971 assert( k >= 0 && k < thedim );
5990 if ( rn3 > thedim*verySparseFactor4right )
5992 ridx3[0] = thedim - 1;
6002 for ( i = j = 0; i < rn3; ++i )
6005 assert( k >= 0 && k < thedim );
6041 Real* rhs2, int* ridx2, int rn2 )
6043 METHOD( "CLUFactor::vSolveRightNoNZ()" );
6045 assert( rn2 >= 0 && rn2 <= thedim );
6049 *ridx2 = thedim - 1;
6061 for ( i = j = 0; i < rn2; ++i )
6064 assert( k >= 0 && k < thedim );
6094 Real* vec, int* idx,
6095 Real* rhs, int* ridx, int rn )
6097 METHOD( "CLUFactor::vSolveLeft()" );
6102 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6106 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6120 Real* vec, int* idx,
6121 Real* rhs, int* ridx, int rn,
6123 Real* rhs2, int* ridx2, int rn2 )
6125 METHOD( "CLUFactor::vSolveLeft2()" );
6130 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6136 rn = solveUleft( eps, vec, idx, rhs, ridx, rn );
6151 Real* rhs2, int* ridx2, int rn2 )
6153 METHOD( "CLUFactor::vSolveLeftNoNZ()" );
|