[Help-glpk] [proposal] spx_update_dvec and spx_update_gvec

[Henri Gourvest]

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;
}