[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-glpk] [proposal] spx_update_dvec and spx_update_gvec
From: |
Henri Gourvest |
Subject: |
[Help-glpk] [proposal] spx_update_dvec and spx_update_gvec |
Date: |
Thu, 16 Jun 2005 20:15:21 +0200 |
User-agent: |
Mozilla Thunderbird 1.0 (Windows/20041206) |
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;
}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Help-glpk] [proposal] spx_update_dvec and spx_update_gvec,
Henri Gourvest <=