qp_solver_impl.h

很多二维 三维几何计算算法 C++ 类库

    }

    // replace variable
    replace_variable( no_ineq);

    // pivot step done
    i = j = -1;
}

template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
replace_variable_original_original( )
{
    // updates for the upper bounded case
    replace_variable_original_original_upd_r(Is_nonnegative());
    
    int  k = in_B[ i];

    // replace original variable [ in: j | out: i ]
    in_B  [ i] = -1;
    in_B  [ j] = k;
       B_O[ k] = j;

    minus_c_B[ k] = 
      ( is_phaseI ? 
	( j < qp_n ? et0 : -aux_c[j-qp_n-slack_A.size()]) : -ET( qp_c[ j]));

    if ( is_phaseI) {
	if ( j >= qp_n) ++art_basic;
	if ( i >= qp_n) --art_basic;
    }

    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) print_basis();
    }
	    
    // update basis inverse
    inv_M_B.enter_original_leave_original( q_x_O.begin(), k);
}

// update of the vector r for U_5 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Standard form      
void  QP_solver<Q, ET, Tags>::
replace_variable_original_original_upd_r(Tag_true )
{
}

// update of the vector r for U_5 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Upper bounded      
void  QP_solver<Q, ET, Tags>::
replace_variable_original_original_upd_r(Tag_false )
{
    ET      x_j, x_i;
    
    if (is_artificial(j)) {
        if (!is_artificial(i)) {
            x_i = (ratio_test_bound_index == LOWER) ? qp_l[i] : qp_u[i];
            update_r_C_r_S_B__i(x_i);
            // update x_O_v_i
            x_O_v_i[i] = ratio_test_bound_index;
        }
    } else {
        x_j = nonbasic_original_variable_value(j);
        if (is_artificial(i)) {
            update_r_C_r_S_B__j(x_j);
        } else {
            x_i = (ratio_test_bound_index == LOWER) ? qp_l[i] : qp_u[i];
            update_r_C_r_S_B__j_i(x_j, x_i);
            // update x_O_v_i
            x_O_v_i[i] = ratio_test_bound_index;
        }
        // update x_O_v_i
        x_O_v_i[j] = BASIC;
    }
}


template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_slack( )
{
    
    // updates for the upper bounded case
    replace_variable_slack_slack_upd_r(Is_nonnegative()); 
    
    int  k = in_B[ i];

    // replace slack variable [ in: j | out: i ]
    in_B  [ i] = -1;
    in_B  [ j] = k;
       B_S[ k] = j;
       S_B[ k] = slack_A[ j-qp_n].first;

    // replace inequality constraint [ in: i | out: j ]
    int old_row = S_B[ k];
    int new_row = slack_A[ i-qp_n].first;
    k = in_C[ old_row];

    in_C[ old_row] = -1;
    in_C[ new_row] = k;
       C[ k      ] = new_row;

     b_C[ k] = ET( qp_b[ new_row]);

    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) print_basis();
    }

    // update basis inverse
    A_row_by_index_accessor  a_accessor( A_accessor( qp_A, 0, qp_n), new_row);
    std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
	       A_row_by_index_iterator( B_O.end  (), a_accessor),
	       tmp_x.begin());
    if ( art_s_i > 0) {                                 // special artificial
	tmp_x[ in_B[ art_s_i]] = ET( art_s[ new_row]);
    }
    inv_M_B.enter_slack_leave_slack( tmp_x.begin(), k);
}

// update of the vector r for U_6 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Standard form      
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_slack_upd_r(Tag_true )
{
}

// update of the vector r for U_6 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Upper bounded      
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_slack_upd_r(Tag_false )
{
    int     sigma_j = slack_A[ j-qp_n].first;
    
    // swap r_gamma_C(sigma_j) in r_C with r_gamma_S_B(sigma_i) in r_S_B
    std::swap(r_C[in_C[sigma_j]], r_S_B[in_B[i]]); 
}


template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_original( )
{
    // updates for the upper bounded case
    replace_variable_slack_original_upd_r(Is_nonnegative()); 
     
    int  k = in_B[ i];

    // leave original variable [ out: i ]
    in_B  [ B_O.back()] = k;
       B_O[ k] = B_O.back();
       in_B  [ i         ] = -1;
       B_O.pop_back();

    minus_c_B[ k] = minus_c_B[ B_O.size()];

    if ( is_phaseI && ( i >= qp_n)) --art_basic;

    // enter slack variable [ in: j ]
    int  old_row = slack_A[ j-qp_n].first;
    in_B  [ j] = B_S.size();
       B_S.push_back( j);
       S_B.push_back( old_row);

    // leave inequality constraint [ out: j ]
    int  l = in_C[ old_row];
     b_C[ l       ] = b_C[ C.size()-1];
       C[ l       ] = C.back();
    in_C[ C.back()] = l;
    in_C[ old_row ] = -1;
       C.pop_back();
    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) print_basis();
    }

    // update basis inverse
    inv_M_B.swap_variable( k);
    inv_M_B.swap_constraint( l);
    inv_M_B.enter_slack_leave_original();
}

// update of the vector r for U_8 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Standard form      
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_original_upd_r(Tag_true )
{
}

// update of the vector r for U_8 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Upper bounded      
void  QP_solver<Q, ET, Tags>::
replace_variable_slack_original_upd_r(Tag_false )
{
    if (!is_artificial(i)) {
        ET  x_i = (ratio_test_bound_index == LOWER) ? qp_l[i] : qp_u[i];
        update_r_C_r_S_B__i(x_i);
    }
    
    int     sigma_j = slack_A[ j-qp_n].first;
    
    // append r_gamma_C(sigma_j) from r_C to r_S_B:
    r_S_B.push_back(r_C[in_C[sigma_j]]);
    
    // remove r_gamma_C(sigma_j) from r_C:
    r_C[in_C[sigma_j]] = r_C.back();
    r_C.pop_back();
    
    // update x_O_v_i
    if (!is_artificial(i)) // original and not artificial?
      x_O_v_i[i] = ratio_test_bound_index;
}


template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
replace_variable_original_slack( )
{
    // updates for the upper bounded case
    replace_variable_original_slack_upd_r(Is_nonnegative());
    
    int  k = in_B[ i];

    // enter original variable [ in: j ]

    minus_c_B[ B_O.size()]
      = ( is_phaseI ? 
	  ( j < qp_n ? et0 : -aux_c[j-qp_n-slack_A.size()]) : -ET( qp_c[ j]));
    

    in_B  [ j] = B_O.size();
       B_O.push_back( j);

    if ( is_phaseI && ( j >= qp_n)) ++art_basic;

    // leave slack variable [ out: i ]
       B_S[ k         ] = B_S.back();
       S_B[ k         ] = S_B.back();
    in_B  [ B_S.back()] = k;
    in_B  [ i         ] = -1; 
       B_S.pop_back();
       S_B.pop_back();

    // enter inequality constraint [ in: i ]
    int new_row = slack_A[ i-qp_n].first;

     b_C[ C.size()] = ET( qp_b[ new_row]);
    in_C[ new_row ] = C.size();
       C.push_back( new_row);
    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) print_basis();
    }

    // update basis inverse
    A_row_by_index_accessor  a_accessor( A_accessor( qp_A, 0, qp_n), new_row);
    std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
	       A_row_by_index_iterator( B_O.end  (), a_accessor),
	       tmp_x.begin());
    if ( art_s_i > 0) {                                 // special art.
	tmp_x[ in_B[ art_s_i]] = ET( art_s[ new_row]);
    }
    inv_M_B.enter_original_leave_slack( q_x_O.begin(), tmp_x.begin());
    
}

// update of the vector r for U_7 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Standard form      
void  QP_solver<Q, ET, Tags>::
replace_variable_original_slack_upd_r(Tag_true )
{
}

// update of the vector r for U_7 with upper bounding, note that we 
// need the headings C, and S_{B} before they are updated
template < typename Q, typename ET, typename Tags >                            // Upper bounded      
void  QP_solver<Q, ET, Tags>::
replace_variable_original_slack_upd_r(Tag_false )
{
    if (!is_artificial(j)) {
        ET x_j = nonbasic_original_variable_value(j);
        update_r_C_r_S_B__j(x_j);
    }
    
    // append r_gamma_S_B(sigma_i) from r_S_B to r_C
    r_C.push_back(r_S_B[in_B[i]]);
    
    // remove r_gamma_S_B(sigma_i) from r_S_B
    r_S_B[in_B[i]] = r_S_B.back();
    r_S_B.pop_back();
    
    // update x_O_v_i
    if (!is_artificial(j)) {
        x_O_v_i[j] = BASIC;
    }
}


template < typename Q, typename ET, typename Tags >
void  QP_solver<Q, ET, Tags>::
remove_artificial_variable_and_constraint( )
{
    // updates for the upper bounded case
    remove_artificial_variable_and_constraint_upd_r(Is_nonnegative());
    
    int  k = in_B[ i];

    // leave artificial (original) variable [ out: i ]
    in_B  [ B_O.back()] = k;
       B_O[ k] = B_O.back();
       in_B  [ i         ] = -1;
       B_O.pop_back();

    minus_c_B[ k] = minus_c_B[ B_O.size()];

    if ( is_phaseI && ( i >= qp_n)) --art_basic;

    int old_row = art_A[i - qp_n - slack_A.size()].first;

    // leave its equality constraint 
    int  l = in_C[ old_row];
     b_C[ l       ] = b_C[ C.size()-1];
       C[ l       ] = C.back();
    in_C[ C.back()] = l;
    in_C[ old_row ] = -1;
       C.pop_back();
    // diagnostic output
    CGAL_qpe_debug {
	if ( vout2.verbose()) print_basis();
    }

    // update basis inverse
    inv_M_B.swap_variable( k);
    inv_M_B.swap_constraint( l);
    inv_M_B.enter_slack_leave_original();
}

// update of the vector r with upper bounding for the removal of an
// artificial variable with its equality constraint, note that we 
// need the headings C before it is updated
template < typename Q, typename ET, typename Tags >                                 // Standard form
void  QP_solver<Q, ET, Tags>::
remove_artificial_variable_and_constraint_upd_r(Tag_true )
{
}

// update of the vector r with upper bounding for the removal of an
// artificial variable with its equality constraint, note that we 
// need the headings C before it is updated
template < typename Q, typename ET, typename Tags >                                 // Upper bounded
void  QP_solver<Q, ET, Tags>::
remove_artificial_variable_and_constraint_upd_r(Tag_false )
{
    int sigma_i = art_A[i - qp_n - slack_A.size()].first;
    
    // remove r_gamma_C(sigma_i) from r_C
    r_C[in_C[sigma_i]] = r_C.back();
    r_C.pop_back();
}

// update that occurs only with upper bounding in ratio test step 1
template < typename Q, typename ET, typename Tags >            
void  QP_solver<Q, ET, Tags>::
enter_and_leave_variable( )
{
    
    CGAL_qpe_assertion((i == j) && (i >= 0));
    
    CGAL_qpe_debug {
	vout2 <<   "<--> nonbasic (" << variable_type( j) << ") variable " << j
	      << " enters and leaves basis" << std::endl << std::endl;
    }

    
    ET diff;
    ET x_j = nonbasic_original_variable_value(j);
    
    if (ratio_test_bound_index == LOWER) {
        diff = x_j - ET(qp_l[j]);
    } else {
        diff = x_j - ET(qp_u[j]);
    }
    
    if (is_phaseI) {
        update_r_C_r_S_B__j(diff);
    } else {
        update_w_r_B_O__j(diff);
        update_r_C_r_S_B__j(diff);
    }
    
    x_O_v_i[j] = ratio_test_bound_index;
    
    // notify pricing strategy (it has called enter_basis on i before)
    strategyP->leaving_basis (i);

    // variable entered and left basis
    i = -1; j = -1;
}



// enter variable into basis
template < typename Q, typename ET, typename Tags >
void
QP_solver<Q, ET, Tags>::
enter_variable( )
{
  CGAL_qpe_assertion (is_phaseII);
    CGAL_qpe_debug {
	vout2 << "--> nonbasic (" << variable_type( j) << ") variable "
	      << j << " enters basis" << std::endl << std::endl;
    }

    // update basis & basis inverse: