自学教程：C++ xfault函数代码示例

void mpl_put_col_value(MPL *mpl, int j, double val){     if (mpl->phase != 3)         xfault("mpl_put_col_value: invalid call sequence/n");      if (!(1 <= j && j <= mpl->n))         xfault(         "mpl_put_col_value: j = %d; column number out of range/n", j);      mpl->col[j]->prim = val;      return;}

int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,      int ind[], double val[]), void *info){     int ret;      if (m < 1)         xfault("fhv_factorize: m = %d; invalid parameter/n", m);      if (m > M_MAX)         xfault("fhv_factorize: m = %d; matrix too big/n", m);      fhv->m = m;      /* invalidate the factorization */      fhv->valid = 0;      /* allocate/reallocate arrays, if necessary */      if (fhv->hh_ind == NULL)         fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int));      if (fhv->hh_ptr == NULL)         fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int));      if (fhv->hh_len == NULL)         fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int));      if (fhv->m_max < m)      {  if (fhv->p0_row != NULL) xfree(fhv->p0_row);         if (fhv->p0_col != NULL) xfree(fhv->p0_col);         if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);         if (fhv->cc_val != NULL) xfree(fhv->cc_val);         fhv->m_max = m + 100;         fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int));         fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int));         fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int));         fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double));      }      /* try to factorize the basis matrix */      switch (luf_factorize(fhv->luf, m, col, info))      {  case 0:            break;         case LUF_ESING:            ret = FHV_ESING;            goto done;         case LUF_ECOND:            ret = FHV_ECOND;            goto done;         default:            xassert(fhv != fhv);      }      /* the basis matrix has been successfully factorized */      fhv->valid = 1;      /* H := I */      fhv->hh_nfs = 0;      /* P0 := P */      memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m);      memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m);      /* currently H has no factors */      fhv->nnz_h = 0;      ret = 0;done: /* return to the calling program */      return ret;}

double mpl_get_row_c0(MPL *mpl, int i){     ELEMCON *con;      double c0;      if (mpl->phase != 3)         xfault("mpl_get_row_c0: invalid call sequence/n");      if (!(1 <= i && i <= mpl->m))         xfault("mpl_get_row_c0: i = %d; row number out of range/n",            i);      con = mpl->row[i];      if (con->con->lbnd == NULL && con->con->ubnd == NULL)         c0 = - con->lbnd;      else         c0 = 0.0;      return c0;}

void lpx_set_class(LPX *lp, int klass){     /* set problem class */      xassert(lp == lp);      if (!(klass == LPX_LP || klass == LPX_MIP))         xfault("lpx_set_class: invalid problem class/n");      return;}

LUX *lux_create(int n){     LUX *lux;      int k;      if (n < 1)         xfault("lux_create: n = %d; invalid parameter/n", n);      lux = xmalloc(sizeof(LUX));      lux->n = n;      lux->pool = dmp_create_poolx(sizeof(LUXELM));      lux->F_row = xcalloc(1+n, sizeof(LUXELM *));      lux->F_col = xcalloc(1+n, sizeof(LUXELM *));      lux->V_piv = xcalloc(1+n, sizeof(mpq_t));      lux->V_row = xcalloc(1+n, sizeof(LUXELM *));      lux->V_col = xcalloc(1+n, sizeof(LUXELM *));      lux->P_row = xcalloc(1+n, sizeof(int));      lux->P_col = xcalloc(1+n, sizeof(int));      lux->Q_row = xcalloc(1+n, sizeof(int));      lux->Q_col = xcalloc(1+n, sizeof(int));      for (k = 1; k <= n; k++)      {  lux->F_row[k] = lux->F_col[k] = NULL;         mpq_init(lux->V_piv[k]);         mpq_set_si(lux->V_piv[k], 1, 1);         lux->V_row[k] = lux->V_col[k] = NULL;         lux->P_row[k] = lux->P_col[k] = k;         lux->Q_row[k] = lux->Q_col[k] = k;      }      lux->rank = n;      return lux;}

/******************************************************************************** Function Name  : HardFault_Handler* Description    : This function handles Hard Fault exception.* Input          : None* Output         : None* Return         : None*******************************************************************************/void HardFault_Handler(void){  xfault();  /* Go to infinite loop when Hard Fault exception occurs */  while (1)  {  }}

int mpl_get_mat_row(MPL *mpl, int i, int ndx[], double val[]){     FORMULA *term;      int len = 0;      if (mpl->phase != 3)         xfault("mpl_get_mat_row: invalid call sequence/n");      if (!(1 <= i && i <= mpl->m))         xfault("mpl_get_mat_row: i = %d; row number out of range/n",            i);      for (term = mpl->row[i]->form; term != NULL; term = term->next)      {  xassert(term->var != NULL);         len++;         xassert(len <= mpl->n);         if (ndx != NULL) ndx[len] = term->var->j;         if (val != NULL) val[len] = term->coef;      }      return len;}

/******************************************************************************** Function Name  : MemManage_Handler* Description    : This function handles Memory Manage exception.* Input          : None* Output         : None* Return         : None*******************************************************************************/void MemManage_Handler(void){  xfault();  /* Go to infinite loop when Memory Manage exception occurs */  while (1)  {  }}

double rng_uniform(RNG *rand, double a, double b){     double x;      if (a >= b)         xfault("rng_uniform: a = %g, b = %g; invalid range/n", a, b);      x = rng_unif_01(rand);      x = a * (1.0 - x) + b * x;      xassert(a <= x && x <= b);      return x;}

void scf_solve_it(SCF *scf, int tr, double x[]){     if (scf->rank < scf->n)         xfault("scf_solve_it: singular matrix/n");      if (!tr)         solve(scf, x);      else         tsolve(scf, x);      return;}

SSX *ssx_create(int m, int n, int nnz){     SSX *ssx;      int i, j, k;      if (m < 1)         xfault("ssx_create: m = %d; invalid number of rows/n", m);      if (n < 1)         xfault("ssx_create: n = %d; invalid number of columns/n", n);      if (nnz < 0)         xfault("ssx_create: nnz = %d; invalid number of non-zero const"            "raint coefficients/n", nnz);      ssx = xmalloc(sizeof(SSX));      ssx->m = m;      ssx->n = n;      ssx->type = xcalloc(1+m+n, sizeof(int));      ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));      for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);      ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));      for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);      ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));      for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);      ssx->A_ptr = xcalloc(1+n+1, sizeof(int));      ssx->A_ptr[n+1] = nnz+1;      ssx->A_ind = xcalloc(1+nnz, sizeof(int));      ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));      for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);      ssx->stat = xcalloc(1+m+n, sizeof(int));      ssx->Q_row = xcalloc(1+m+n, sizeof(int));      ssx->Q_col = xcalloc(1+m+n, sizeof(int));      ssx->binv = bfx_create_binv();      ssx->bbar = xcalloc(1+m, sizeof(mpq_t));      for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);      ssx->pi = xcalloc(1+m, sizeof(mpq_t));      for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);      ssx->cbar = xcalloc(1+n, sizeof(mpq_t));      for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);      ssx->rho = xcalloc(1+m, sizeof(mpq_t));      for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);      ssx->ap = xcalloc(1+n, sizeof(mpq_t));      for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);      ssx->aq = xcalloc(1+m, sizeof(mpq_t));      for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);      mpq_init(ssx->delta);      return ssx;}

int mpl_read_model(MPL *mpl, char *file, int skip_data){     if (mpl->phase != 0)         xfault("mpl_read_model: invalid call sequence/n");      if (file == NULL)         xfault("mpl_read_model: no input filename specified/n");      /* set up error handler */      if (setjmp(mpl->jump)) goto done;      /* translate model section */      mpl->phase = 1;      xprintf("Reading model section from %s.../n", file);      open_input(mpl, file);      model_section(mpl);      if (mpl->model == NULL)         error(mpl, "empty model section not allowed");      /* save name of the input text file containing model section for         error diagnostics during the generation phase */      mpl->mod_file = xcalloc(strlen(file)+1, sizeof(char));      strcpy(mpl->mod_file, mpl->in_file);      /* allocate content arrays for all model objects */      alloc_content(mpl);      /* optional data section may begin with the keyword 'data' */      if (is_keyword(mpl, "data"))      {  if (skip_data)         {  warning(mpl, "data section ignored");            goto skip;         }         mpl->flag_d = 1;         get_token(mpl /* data */);         if (mpl->token != T_SEMICOLON)            error(mpl, "semicolon missing where expected");         get_token(mpl /* ; */);         /* translate data section */         mpl->phase = 2;         xprintf("Reading data section from %s.../n", file);         data_section(mpl);      }      /* process end statement */      end_statement(mpl);skip: xprintf("%d line%s were read/n",         mpl->line, mpl->line == 1 ? "" : "s");      close_input(mpl);done: /* return to the calling program */      return mpl->phase;}

int mpl_get_col_kind(MPL *mpl, int j){     int kind;      if (mpl->phase != 3)         xfault("mpl_get_col_kind: invalid call sequence/n");      if (!(1 <= j && j <= mpl->n))         xfault("mpl_get_col_kind: j = %d; column number out of range/n"            , j);      switch (mpl->col[j]->var->type)      {  case A_NUMERIC:            kind = MPL_NUM; break;         case A_INTEGER:            kind = MPL_INT; break;         case A_BINARY:            kind = MPL_BIN; break;         default:            xassert(mpl != mpl);      }      return kind;}

int mpl_get_row_kind(MPL *mpl, int i){     int kind;      if (mpl->phase != 3)         xfault("mpl_get_row_kind: invalid call sequence/n");      if (!(1 <= i && i <= mpl->m))         xfault("mpl_get_row_kind: i = %d; row number out of range/n",            i);      switch (mpl->row[i]->con->type)      {  case A_CONSTRAINT:            kind = MPL_ST; break;         case A_MINIMIZE:            kind = MPL_MIN; break;         case A_MAXIMIZE:            kind = MPL_MAX; break;         default:            xassert(mpl != mpl);      }      return kind;}

int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],      const double val[]){     int ret;      if (!bfd->valid)         xfault("bfd_update_it: the factorization is not valid/n");      /* try to update the factorization */      if (bfd->fhv != NULL)      {  switch (fhv_update_it(bfd->fhv, j, len, ind, val))         {  case 0:               break;            case FHV_ESING:               bfd->valid = 0;               ret = BFD_ESING;               goto done;            case FHV_ECHECK:               bfd->valid = 0;               ret = BFD_ECHECK;               goto done;            case FHV_ELIMIT:               bfd->valid = 0;               ret = BFD_ELIMIT;               goto done;            case FHV_EROOM:               bfd->valid = 0;               ret = BFD_EROOM;               goto done;            default:               xassert(bfd != bfd);         }      }      else if (bfd->lpf != NULL)      {  switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val))         {  case 0:               break;            case LPF_ESING:               bfd->valid = 0;               ret = BFD_ESING;               goto done;            case LPF_ELIMIT:               bfd->valid = 0;               ret = BFD_ELIMIT;               goto done;            default:               xassert(bfd != bfd);         }      }      else         xassert(bfd != bfd);      /* the factorization has been successfully updated */      /* increase the update count */      bfd->upd_cnt++;      ret = 0;done: /* return to the calling program */      return ret;}

void lpf_btran(LPF *lpf, double x[]){     int m0 = lpf->m0;      int m = lpf->m;      int n = lpf->n;      int *P_row = lpf->P_row;      int *Q_row = lpf->Q_row;      double *fg = lpf->work1;      double *f = fg;      double *g = fg + m0;      int i, ii;#if GLPLPF_DEBUG      double *b;#endif      if (!lpf->valid)         xfault("lpf_btran: the factorization is not valid/n");      xassert(0 <= m && m <= m0 + n);#if GLPLPF_DEBUG      /* save the right-hand side vector */      b = xcalloc(1+m, sizeof(double));      for (i = 1; i <= m; i++) b[i] = x[i];#endif      /* (f g) := Q * (b 0) */      for (i = 1; i <= m0 + n; i++)         fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0);      /* f1 := inv(U'0) * f */#if 0 /* 06/VI-2013 */      luf_v_solve(lpf->luf, 1, f);#else      {  double *work = lpf->lufint->sgf->work;         luf_vt_solve(lpf->lufint->luf, f, work);         memcpy(&f[1], &work[1], m0 * sizeof(double));      }#endif      /* g1 := inv(C') * (g - R' * f1) */      rt_prod(lpf, g, -1.0, f);      scf_solve_it(lpf->scf, 1, g);      /* g2 := g1 */      g = g;      /* f2 := inv(L'0) * (f1 - S' * g2) */      st_prod(lpf, f, -1.0, g);#if 0 /* 06/VI-2013 */      luf_f_solve(lpf->luf, 1, f);#else      luf_ft_solve(lpf->lufint->luf, f);#endif      /* (x y) := P * (f2 g2) */      for (i = 1; i <= m; i++)         x[i] = fg[P_row[i]];#if GLPLPF_DEBUG      /* check relative error in solution */      check_error(lpf, 1, x, b);      xfree(b);#endif      return;}

void bfd_btran(BFD *bfd, double x[]){     if (!bfd->valid)         xfault("bfd_btran: the factorization is not valid/n");      if (bfd->fhv != NULL)         fhv_btran(bfd->fhv, x);      else if (bfd->lpf != NULL)         lpf_btran(bfd->lpf, x);      else         xassert(bfd != bfd);      return;}

char *mpl_get_row_name(MPL *mpl, int i){     char *name = mpl->mpl_buf, *t;      int len;      if (mpl->phase != 3)         xfault("mpl_get_row_name: invalid call sequence/n");      if (!(1 <= i && i <= mpl->m))         xfault("mpl_get_row_name: i = %d; row number out of range/n",            i);      strcpy(name, mpl->row[i]->con->name);      len = strlen(name);      xassert(len <= 255);      t = format_tuple(mpl, '[', mpl->row[i]->memb->tuple);      while (*t)      {  if (len == 255) break;         name[len++] = *t++;      }      name[len] = '/0';      if (len == 255) strcpy(name+252, "...");      xassert(strlen(name) <= 255);      return name;}

int mpl_get_row_bnds(MPL *mpl, int i, double *_lb, double *_ub){     ELEMCON *con;      int type;      double lb, ub;      if (mpl->phase != 3)         xfault("mpl_get_row_bnds: invalid call sequence/n");      if (!(1 <= i && i <= mpl->m))         xfault("mpl_get_row_bnds: i = %d; row number out of range/n",            i);      con = mpl->row[i];#if 0 /* 21/VII-2006 */      if (con->con->lbnd == NULL && con->con->ubnd == NULL)         type = MPL_FR, lb = ub = 0.0;      else if (con->con->ubnd == NULL)         type = MPL_LO, lb = con->lbnd, ub = 0.0;      else if (con->con->lbnd == NULL)         type = MPL_UP, lb = 0.0, ub = con->ubnd;      else if (con->con->lbnd != con->con->ubnd)         type = MPL_DB, lb = con->lbnd, ub = con->ubnd;      else         type = MPL_FX, lb = ub = con->lbnd;#else      lb = (con->con->lbnd == NULL ? -DBL_MAX : con->lbnd);      ub = (con->con->ubnd == NULL ? +DBL_MAX : con->ubnd);      if (lb == -DBL_MAX && ub == +DBL_MAX)         type = MPL_FR, lb = ub = 0.0;      else if (ub == +DBL_MAX)         type = MPL_LO, ub = 0.0;      else if (lb == -DBL_MAX)         type = MPL_UP, lb = 0.0;      else if (con->con->lbnd != con->con->ubnd)         type = MPL_DB;      else         type = MPL_FX;#endif      if (_lb != NULL) *_lb = lb;      if (_ub != NULL) *_ub = ub;      return type;}

int lpx_transform_col(LPX *lp, int len, int ind[], double val[]){     int i, m, t;      double *a, *alfa;      if (!lpx_is_b_avail(lp))         xfault("lpx_transform_col: LP basis is not available/n");      m = lpx_get_num_rows(lp);      /* unpack the column to be transformed to the array a */      a = xcalloc(1+m, sizeof(double));      for (i = 1; i <= m; i++) a[i] = 0.0;      if (!(0 <= len && len <= m))         xfault("lpx_transform_col: len = %d; invalid column length/n",            len);      for (t = 1; t <= len; t++)      {  i = ind[t];         if (!(1 <= i && i <= m))            xfault("lpx_transform_col: ind[%d] = %d; row index out of r"               "ange/n", t, i);         if (val[t] == 0.0)            xfault("lpx_transform_col: val[%d] = 0; zero coefficient no"               "t allowed/n", t);         if (a[i] != 0.0)            xfault("lpx_transform_col: ind[%d] = %d; duplicate row indi"               "ces not allowed/n", t, i);         a[i] = val[t];      }      /* solve the system B*a = alfa to compute the vector alfa */      alfa = a, lpx_ftran(lp, alfa);      /* store resultant coefficients */      len = 0;      for (i = 1; i <= m; i++)      {  if (alfa[i] != 0.0)         {  len++;            ind[len] = lpx_get_b_info(lp, i);            val[len] = alfa[i];         }      }      xfree(a);      return len;}

int mpl_postsolve(MPL *mpl){     if (!(mpl->phase == 3 && !mpl->flag_p))         xfault("mpl_postsolve: invalid call sequence/n");      /* set up error handler */      if (setjmp(mpl->jump)) goto done;      /* perform postsolving */      postsolve_model(mpl);      flush_output(mpl);      /* postsolving phase has been finished */      xprintf("Model has been successfully processed/n");done: /* return to the calling program */      return mpl->phase;}

int mpl_get_col_bnds(MPL *mpl, int j, double *_lb, double *_ub){     ELEMVAR *var;      int type;      double lb, ub;      if (mpl->phase != 3)         xfault("mpl_get_col_bnds: invalid call sequence/n");      if (!(1 <= j && j <= mpl->n))         xfault("mpl_get_col_bnds: j = %d; column number out of range/n"            , j);      var = mpl->col[j];#if 0 /* 21/VII-2006 */      if (var->var->lbnd == NULL && var->var->ubnd == NULL)         type = MPL_FR, lb = ub = 0.0;      else if (var->var->ubnd == NULL)         type = MPL_LO, lb = var->lbnd, ub = 0.0;      else if (var->var->lbnd == NULL)         type = MPL_UP, lb = 0.0, ub = var->ubnd;      else if (var->var->lbnd != var->var->ubnd)         type = MPL_DB, lb = var->lbnd, ub = var->ubnd;      else         type = MPL_FX, lb = ub = var->lbnd;#else      lb = (var->var->lbnd == NULL ? -DBL_MAX : var->lbnd);      ub = (var->var->ubnd == NULL ? +DBL_MAX : var->ubnd);      if (lb == -DBL_MAX && ub == +DBL_MAX)         type = MPL_FR, lb = ub = 0.0;      else if (ub == +DBL_MAX)         type = MPL_LO, ub = 0.0;      else if (lb == -DBL_MAX)         type = MPL_UP, lb = 0.0;      else if (var->var->lbnd != var->var->ubnd)         type = MPL_DB;      else         type = MPL_FX;#endif      if (_lb != NULL) *_lb = lb;      if (_ub != NULL) *_ub = ub;      return type;}

void lux_solve(LUX *lux, int tr, mpq_t x[]){     if (lux->rank < lux->n)         xfault("lux_solve: LU-factorization has incomplete rank/n");      if (!tr)      {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */         lux_f_solve(lux, 0, x);         lux_v_solve(lux, 0, x);      }      else      {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */         lux_v_solve(lux, 1, x);         lux_f_solve(lux, 1, x);      }      return;}

void lpx_eval_b_dual(LPX *lp, double row_dual[], double col_dual[]){     int i, j, k, m, n, len, *ind;      double dj, *cB, *pi, *val;      if (!lpx_is_b_avail(lp))         xfault("lpx_eval_b_dual: LP basis is not available/n");      m = lpx_get_num_rows(lp);      n = lpx_get_num_cols(lp);      /* store zero reduced costs of basic auxiliary and structural         variables and build the vector cB of objective coefficients at         basic variables */      cB = xcalloc(1+m, sizeof(double));      for (i = 1; i <= m; i++)      {  k = lpx_get_b_info(lp, i);         /* xB[i] is k-th original variable */         xassert(1 <= k && k <= m+n);         if (k <= m)         {  row_dual[k] = 0.0;            cB[i] = 0.0;         }         else         {  col_dual[k-m] = 0.0;            cB[i] = lpx_get_obj_coef(lp, k-m);         }      }      /* solve the system B'*pi = cB to compute the vector pi */      pi = cB, lpx_btran(lp, pi);      /* compute reduced costs of non-basic auxiliary variables */      for (i = 1; i <= m; i++)      {  if (lpx_get_row_stat(lp, i) != LPX_BS)            row_dual[i] = - pi[i];      }      /* compute reduced costs of non-basic structural variables */      ind = xcalloc(1+m, sizeof(int));      val = xcalloc(1+m, sizeof(double));      for (j = 1; j <= n; j++)      {  if (lpx_get_col_stat(lp, j) != LPX_BS)         {  dj = lpx_get_obj_coef(lp, j);            len = lpx_get_mat_col(lp, j, ind, val);            for (k = 1; k <= len; k++) dj += val[k] * pi[ind[k]];            col_dual[j] = dj;         }      }      xfree(ind);      xfree(val);      xfree(cB);      return;}

void fhv_btran(FHV *fhv, double x[]){     int *pp_row = fhv->luf->pp_row;      int *pp_col = fhv->luf->pp_col;      int *p0_row = fhv->p0_row;      int *p0_col = fhv->p0_col;      if (!fhv->valid)         xfault("fhv_btran: the factorization is not valid/n");      /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */      luf_v_solve(fhv->luf, 1, x);      fhv_h_solve(fhv, 1, x);      fhv->luf->pp_row = p0_row;      fhv->luf->pp_col = p0_col;      luf_f_solve(fhv->luf, 1, x);      fhv->luf->pp_row = pp_row;      fhv->luf->pp_col = pp_col;      return;}

void lpf_ftran(LPF *lpf, double x[]){     int m0 = lpf->m0;      int m = lpf->m;      int n  = lpf->n;      int *P_col = lpf->P_col;      int *Q_col = lpf->Q_col;      double *fg = lpf->work1;      double *f = fg;      double *g = fg + m0;      int i, ii;#if _GLPLPF_DEBUG      double *b;#endif      if (!lpf->valid)         xfault("lpf_ftran: the factorization is not valid/n");      xassert(0 <= m && m <= m0 + n);#if _GLPLPF_DEBUG      /* save the right-hand side vector */      b = xcalloc(1+m, sizeof(double));      for (i = 1; i <= m; i++) b[i] = x[i];#endif      /* (f g) := inv(P) * (b 0) */      for (i = 1; i <= m0 + n; i++)         fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0);      /* f1 := inv(L0) * f */      luf_f_solve(lpf->luf, 0, f);      /* g1 := g - S * f1 */      s_prod(lpf, g, -1.0, f);      /* g2 := inv(C) * g1 */      scf_solve_it(lpf->scf, 0, g);      /* f2 := inv(U0) * (f1 - R * g2) */      r_prod(lpf, f, -1.0, g);      luf_v_solve(lpf->luf, 0, f);      /* (x y) := inv(Q) * (f2 g2) */      for (i = 1; i <= m; i++)         x[i] = fg[Q_col[i]];#if _GLPLPF_DEBUG      /* check relative error in solution */      check_error(lpf, 0, x, b);      xfree(b);#endif      return;}