Hi all,
"spx_update_dvec" and "spx_update_gvec" are precision sensitives
regarding to speed and quality results.
I've modified these methods and It speed up most of my problems, for
example "jssp.mod" is solved in 38sec instead of 61 sec.
Any feedback ?
Henri Gourvest
http://www.progdigy.com
void spx_update_dvec(SPX *spx)
{ int m = spx->m;
int n = spx->n;
int *typx = spx->typx;
int *AT_ptr = spx->AT_ptr;
int *AT_ind = spx->AT_ind;
double *AT_val = spx->AT_val;
int *indx = spx->indx;
int p = spx->p;
int q = spx->q;
double *ap = spx->ap;
double *aq = spx->aq;
double *dvec = spx->dvec;
int *refsp = spx->refsp;
double *w = spx->work;
int i, j, j_beg, j_end, j_ptr, k, ref_k, ref_p, ref_q;
double ap_j, aq_p, aq_i;
double s_i, t1, sum, temp;
insist(1 <= p && p <= m);
insist(1 <= q && q <= n);
/* check if it's time to reset the reference space */
if (spx->count <= 0)
{ /* yes, do that */
spx_reset_refsp(spx);
goto done;
}
else
{ /* otherwise decrease the count */
spx->count--;
}
/* compute t1 */
t1 = 0.0;
for (j = 1; j <= n; j++)
{ if (j != q && refsp[indx[m+j]])
{ ap_j = ap[j];
t1 += ap_j * ap_j;
}
}
/* compute w */
for (i = 1; i <= m; i++) w[i] = 0.0;
for (j = 1; j <= n; j++)
{ if (j != q && refsp[indx[m+j]])
{ /* w += ap[j] * N[j] */
ap_j = ap[j];
if (ap_j == 0.0) continue;
k = indx[m+j]; /* x[k] = xN[j] */
if (k <= m)
/* xN[j] is auxiliary variable */
w[k] += ap_j;
else
{ /* xN[j] is structural variable */
j_beg = AT_ptr[k-m], j_end = AT_ptr[k-m+1];
for (j_ptr = j_beg; j_ptr < j_end; j_ptr++)
w[AT_ind[j_ptr]] -= ap_j * AT_val[j_ptr];
}
}
}
spx_ftran(spx, w, 0);
/* update the vector delta */
ref_p = refsp[indx[p]]; /* if xB[p] belongs to the space */
ref_q = refsp[indx[m+q]]; /* if xN[q] belongs to the space */
aq_p = aq[p];
insist(aq_p != 0.0);
for (i = 1; i <= m; i++)
{ /* dvec[p] will be computed later */
if (i == p) continue;
/* if xB[i] is free variable, its weight is not used, because
free variables never leave the basis */
k = indx[i];
if (typx[k] == LPX_FR)
{ dvec[i] = 1.0;
continue;
}
if (aq[i] != 0.0){
s_i = dvec[i];
aq_i = aq[i];
if (ref_q != 0.0)
s_i = s_i - (aq_i * aq_i);
if (ref_p != 0.0)
s_i = s_i + aq_i*((2.0 * w[i]) + (t1*aq_i) + aq_i)/(aq_p *
aq_p);
else
s_i = s_i + aq_i*((2.0 * w[i]) + (t1*aq_i))/(aq_p * aq_p);
/* update dvec[i] */
if (s_i < DBL_EPSILON)
dvec[i] = 1.0; else
dvec[i] = s_i;
}
}
/* compute exact value of dvec[p] */
sum = (ref_q ? 1.0 : 0.0);
double temp2 = aq_p * aq_p;
temp = 0.0;
for(j = 1; j <= n; j++){
if (j == q){
if (ref_p != 0)
temp += 1.0;
}
else
if (refsp[indx[m + j]] != 0){
ap_j = ap[j];
temp += (ap_j * ap_j);
}
}
dvec[p] = sum + temp / temp2;
done: return;
}
void spx_update_gvec(SPX *spx)
{ int m = spx->m;
int n = spx->n;
int *AT_ptr = spx->AT_ptr;
int *AT_ind = spx->AT_ind;
double *AT_val = spx->AT_val;
int *tagx = spx->tagx;
int *indx = spx->indx;
int p = spx->p;
int q = spx->q;
double *ap = spx->ap;
double *aq = spx->aq;
double *gvec = spx->gvec;
int *refsp = spx->refsp;
int i, j, j_beg, j_end, j_ptr, k, ref_k, ref_p, ref_q;
double ap_j, ap_q, aq_i;
double s_j, t1, t2, t3, sum, temp;
double *w = spx->work;
insist(1 <= p && p <= m);
insist(1 <= q && q <= n);
/* check if it's time to reset the reference space */
if (spx->count <= 0)
{ /* yes, do that */
spx_reset_refsp(spx);
goto done;
}
else
{ /* otherwise decrease the count */
spx->count--;
}
/* compute t1 and w */
t1 = 0.0;
for (i = 1; i <= m; i++)
{ if (i != p && refsp[indx[i]])
{ w[i] = aq[i];
t1 += w[i] * w[i];
}
else
w[i] = 0.0;
}
spx_btran(spx, w);
/* update the vector gamma */
ref_p = refsp[indx[p]]; /* if xB[p] belongs to the space */
ref_q = refsp[indx[m+q]]; /* if xN[q] belongs to the space */
ap_q = ap[q];
insist(ap_q != 0.0);
for (j = 1; j <= n; j++)
{ /* gvec[q] will be computed later */
if (j == q) continue;
/* if xN[j] is fixed variable, its weight is not used, because
fixed variables never enter the basis */
k = indx[m+j]; /* x[k] = xN[j] */
if (tagx[k] == LPX_NS)
{ gvec[j] = 1.0;
continue;
}
ref_k = refsp[k]; /* if xN[j] belongs to the space */
/* compute s[j] */
ap_j = ap[j];
s_j = gvec[j];
if (ap[j] != 0.0){
ap_j = ap[j];
if (ref_p != 0)
s_j = s_j - (ap_j * ap_j);
/* t2 := N[j] * w */
if (k <= m){
/* x[k] is auxiliary variable */
t2 = +w[k];
}
else
{
/* x[k] is structural variable */
t2 = 0.0;
j_beg = AT_ptr[k - m];
j_end = AT_ptr[k - m + 1];
j_ptr = j_beg;
while (j_ptr < j_end){
t2 = t2 - (AT_val[j_ptr] * w[AT_ind[j_ptr]]);
j_ptr++;
}
}
if (ref_q != 0)
s_j = s_j + ap_j * ((2 * t2 * ap_q) + (t1 * ap_j) + ap_j) /
(ap_q * ap_q);
else
s_j = s_j + ap_j * ((2 * t2 * ap_q) + (t1 * ap_j)) / (ap_q
* ap_q);
/* update gvec[j] */
if (s_j < DBL_EPSILON)
gvec[j] = 1.0;
else
gvec[j] = s_j;
}
}
/* compute exact value of gvec[q] */
sum = (ref_p ? 1.0 : 0.0);
temp = 0.0;
double temp2 = ap_q * ap_q;
for(i = 1; i <= m; i++){
if (i == p){
if (ref_q != 0)
temp += 1.0;
}
else
{
if (refsp[indx[i]] != 0){
aq_i = aq[i];
temp = temp + (aq_i * aq_i);
}
}
}
gvec[q] = sum + temp / temp2;
done: return;
}
_______________________________________________
Help-glpk mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-glpk
