glplpx04.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

            lpx_get_col_info(lp, j-m, NULL, &xS_j, NULL);
            xS_j /= lpx_get_sjj(lp, j-m);
#endif
            /* g[i] := g[i] - a[i,j] * xS[j] */
            rii = lpx_get_rii(lp, i);
            sjj = lpx_get_sjj(lp, j-m);
            g_i -= (rii * A_val[t] * sjj) * xS_j;
         }
         /* unscale xR[i] and g[i] (if required) */
         if (!scaled)
         {  rii = lpx_get_rii(lp, i);
            xR_i /= rii, g_i /= rii;
         }
         /* determine absolute error */
         temp = fabs(g_i);
         if (kkt->pe_ae_max < temp)
            kkt->pe_ae_max = temp, kkt->pe_ae_row = i;
         /* determine relative error */
         temp /= (1.0 + fabs(xR_i));
         if (kkt->pe_re_max < temp)
            kkt->pe_re_max = temp, kkt->pe_re_row = i;
      }
      xfree(A_ndx);
      xfree(A_val);
      /* estimate the solution quality */
      if (kkt->pe_re_max <= 1e-9)
         kkt->pe_quality = 'H';
      else if (kkt->pe_re_max <= 1e-6)
         kkt->pe_quality = 'M';
      else if (kkt->pe_re_max <= 1e-3)
         kkt->pe_quality = 'L';
      else
         kkt->pe_quality = '?';
      /*--------------------------------------------------------------*/
      /* compute largest absolute and relative errors and corresponding
         variable indices for the condition (KKT.PB) */
      kkt->pb_ae_max = 0.0, kkt->pb_ae_ind = 0;
      kkt->pb_re_max = 0.0, kkt->pb_re_ind = 0;
      for (k = 1; k <= m+n; k++)
      {  /* determine x[k] */
         if (k <= m)
         {  lpx_get_row_bnds(lp, k, &typx, &lb, &ub);
            rii = lpx_get_rii(lp, k);
            lb *= rii;
            ub *= rii;
            lpx_get_row_info(lp, k, &tagx, &x_k, NULL);
            x_k *= rii;
         }
         else
         {  lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub);
            sjj = lpx_get_sjj(lp, k-m);
            lb /= sjj;
            ub /= sjj;
            lpx_get_col_info(lp, k-m, &tagx, &x_k, NULL);
            x_k /= sjj;
         }
         /* skip non-basic variable */
         if (tagx != LPX_BS) continue;
         /* compute h[k] */
         h_k = 0.0;
         switch (typx)
         {  case LPX_FR:
               break;
            case LPX_LO:
               if (x_k < lb) h_k = x_k - lb;
               break;
            case LPX_UP:
               if (x_k > ub) h_k = x_k - ub;
               break;
            case LPX_DB:
            case LPX_FX:
               if (x_k < lb) h_k = x_k - lb;
               if (x_k > ub) h_k = x_k - ub;
               break;
            default:
               xassert(typx != typx);
         }
         /* unscale x[k] and h[k] (if required) */
         if (!scaled)
         {  if (k <= m)
            {  rii = lpx_get_rii(lp, k);
               x_k /= rii, h_k /= rii;
            }
            else
            {  sjj = lpx_get_sjj(lp, k-m);
               x_k *= sjj, h_k *= sjj;
            }
         }
         /* determine absolute error */
         temp = fabs(h_k);
         if (kkt->pb_ae_max < temp)
            kkt->pb_ae_max = temp, kkt->pb_ae_ind = k;
         /* determine relative error */
         temp /= (1.0 + fabs(x_k));
         if (kkt->pb_re_max < temp)
            kkt->pb_re_max = temp, kkt->pb_re_ind = k;
      }
      /* estimate the solution quality */
      if (kkt->pb_re_max <= 1e-9)
         kkt->pb_quality = 'H';
      else if (kkt->pb_re_max <= 1e-6)
         kkt->pb_quality = 'M';
      else if (kkt->pb_re_max <= 1e-3)
         kkt->pb_quality = 'L';
      else
         kkt->pb_quality = '?';
      /*--------------------------------------------------------------*/
      /* compute largest absolute and relative errors and corresponding
         column indices for the condition (KKT.DE) */
      kkt->de_ae_max = 0.0, kkt->de_ae_col = 0;
      kkt->de_re_max = 0.0, kkt->de_re_col = 0;
      A_ndx = xcalloc(1+m, sizeof(int));
      A_val = xcalloc(1+m, sizeof(double));
      for (j = m+1; j <= m+n; j++)
      {  /* determine cS[j] */
#if 0 /* 21/XII-2003 */
         cS_j = coef[j];
#else
         sjj = lpx_get_sjj(lp, j-m);
         cS_j = lpx_get_obj_coef(lp, j-m) * sjj;
#endif
         /* determine dS[j] */
#if 0 /* 21/XII-2003 */
         if (tagx[j] == LPX_BS)
            dS_j = 0.0;
         else
            dS_j = cbar[posx[j] - m];
#else
         lpx_get_col_info(lp, j-m, NULL, NULL, &dS_j);
         dS_j *= sjj;
#endif
         /* u[j] := dS[j] - cS[j] */
         u_j = dS_j - cS_j;
         /* u[j] := u[j] + (j-th column of A) * (dR - cR) */
         beg = 1;
         end = lpx_get_mat_col(lp, j-m, A_ndx, A_val);
         for (t = beg; t <= end; t++)
         {  i = A_ndx[t]; /* a[i,j] != 0 */
            /* determine cR[i] */
#if 0 /* 21/XII-2003 */
            cR_i = coef[i];
#else
            cR_i = 0.0;
#endif
            /* determine dR[i] */
#if 0 /* 21/XII-2003 */
            if (tagx[i] == LPX_BS)
               dR_i = 0.0;
            else
               dR_i = cbar[posx[i] - m];
#else
            lpx_get_row_info(lp, i, NULL, NULL, &dR_i);
            rii = lpx_get_rii(lp, i);
            dR_i /= rii;
#endif
            /* u[j] := u[j] + a[i,j] * (dR[i] - cR[i]) */
            rii = lpx_get_rii(lp, i);
            sjj = lpx_get_sjj(lp, j-m);
            u_j += (rii * A_val[t] * sjj) * (dR_i - cR_i);
         }
         /* unscale cS[j], dS[j], and u[j] (if required) */
         if (!scaled)
         {  sjj = lpx_get_sjj(lp, j-m);
            cS_j /= sjj, dS_j /= sjj, u_j /= sjj;
         }
         /* determine absolute error */
         temp = fabs(u_j);
         if (kkt->de_ae_max < temp)
            kkt->de_ae_max = temp, kkt->de_ae_col = j - m;
         /* determine relative error */
         temp /= (1.0 + fabs(dS_j - cS_j));
         if (kkt->de_re_max < temp)
            kkt->de_re_max = temp, kkt->de_re_col = j - m;
      }
      xfree(A_ndx);
      xfree(A_val);
      /* estimate the solution quality */
      if (kkt->de_re_max <= 1e-9)
         kkt->de_quality = 'H';
      else if (kkt->de_re_max <= 1e-6)
         kkt->de_quality = 'M';
      else if (kkt->de_re_max <= 1e-3)
         kkt->de_quality = 'L';
      else
         kkt->de_quality = '?';
      /*--------------------------------------------------------------*/
      /* compute largest absolute and relative errors and corresponding
         variable indices for the condition (KKT.DB) */
      kkt->db_ae_max = 0.0, kkt->db_ae_ind = 0;
      kkt->db_re_max = 0.0, kkt->db_re_ind = 0;
      for (k = 1; k <= m+n; k++)
      {  /* determine c[k] */
#if 0 /* 21/XII-2003 */
         c_k = coef[k];
#else
         if (k <= m)
            c_k = 0.0;
         else
         {  sjj = lpx_get_sjj(lp, k-m);
            c_k = lpx_get_obj_coef(lp, k-m) / sjj;
         }
#endif
         /* determine d[k] */
#if 0 /* 21/XII-2003 */
         d_k = cbar[j-m];
#else
         if (k <= m)
         {  lpx_get_row_info(lp, k, &tagx, NULL, &d_k);
            rii = lpx_get_rii(lp, k);
            d_k /= rii;
         }
         else
         {  lpx_get_col_info(lp, k-m, &tagx, NULL, &d_k);
            sjj = lpx_get_sjj(lp, k-m);
            d_k *= sjj;
         }
#endif
         /* skip basic variable */
         if (tagx == LPX_BS) continue;
         /* compute v[k] */
         v_k = 0.0;
         switch (tagx)
         {  case LPX_NL:
               switch (dir)
               {  case LPX_MIN:
                     if (d_k < 0.0) v_k = d_k;
                     break;
                  case LPX_MAX:
                     if (d_k > 0.0) v_k = d_k;
                     break;
                  default:
                     xassert(dir != dir);
               }
               break;
            case LPX_NU:
               switch (dir)
               {  case LPX_MIN:
                     if (d_k > 0.0) v_k = d_k;
                     break;
                  case LPX_MAX:
                     if (d_k < 0.0) v_k = d_k;
                     break;
                  default:
                     xassert(dir != dir);
               }
               break;
            case LPX_NF:
               v_k = d_k;
               break;
            case LPX_NS:
               break;
            default:
               xassert(tagx != tagx);
         }
         /* unscale c[k], d[k], and v[k] (if required) */
         if (!scaled)
         {  if (k <= m)
            {  rii = lpx_get_rii(lp, k);
               c_k *= rii, d_k *= rii, v_k *= rii;
            }
            else
            {  sjj = lpx_get_sjj(lp, k-m);
               c_k /= sjj, d_k /= sjj, v_k /= sjj;
            }
         }
         /* determine absolute error */
         temp = fabs(v_k);
         if (kkt->db_ae_max < temp)
            kkt->db_ae_max = temp, kkt->db_ae_ind = k;
         /* determine relative error */
         temp /= (1.0 + fabs(d_k - c_k));
         if (kkt->db_re_max < temp)
            kkt->db_re_max = temp, kkt->db_re_ind = k;
      }
      /* estimate the solution quality */
      if (kkt->db_re_max <= 1e-9)
         kkt->db_quality = 'H';
      else if (kkt->db_re_max <= 1e-6)
         kkt->db_quality = 'M';
      else if (kkt->db_re_max <= 1e-3)
         kkt->db_quality = 'L';
      else
         kkt->db_quality = '?';
      /* complementary slackness is always satisfied by definition for
         any basic solution, so not checked */
      kkt->cs_ae_max = 0.0, kkt->cs_ae_ind = 0;
      kkt->cs_re_max = 0.0, kkt->cs_re_ind = 0;
      kkt->cs_quality = 'H';
      return;
}

/* eof */