1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glplpf.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,978 @@
     1.4 +/* glplpf.c (LP basis factorization, Schur complement version) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glplpf.h"
    1.29 +#include "glpenv.h"
    1.30 +#define xfault xerror
    1.31 +
    1.32 +#define _GLPLPF_DEBUG 0
    1.33 +
    1.34 +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
    1.35 +
    1.36 +#define M_MAX 100000000 /* = 100*10^6 */
    1.37 +/* maximal order of the basis matrix */
    1.38 +
    1.39 +/***********************************************************************
    1.40 +*  NAME
    1.41 +*
    1.42 +*  lpf_create_it - create LP basis factorization
    1.43 +*
    1.44 +*  SYNOPSIS
    1.45 +*
    1.46 +*  #include "glplpf.h"
    1.47 +*  LPF *lpf_create_it(void);
    1.48 +*
    1.49 +*  DESCRIPTION
    1.50 +*
    1.51 +*  The routine lpf_create_it creates a program object, which represents
    1.52 +*  a factorization of LP basis.
    1.53 +*
    1.54 +*  RETURNS
    1.55 +*
    1.56 +*  The routine lpf_create_it returns a pointer to the object created. */
    1.57 +
    1.58 +LPF *lpf_create_it(void)
    1.59 +{     LPF *lpf;
    1.60 +#if _GLPLPF_DEBUG
    1.61 +      xprintf("lpf_create_it: warning: debug mode enabled\n");
    1.62 +#endif
    1.63 +      lpf = xmalloc(sizeof(LPF));
    1.64 +      lpf->valid = 0;
    1.65 +      lpf->m0_max = lpf->m0 = 0;
    1.66 +      lpf->luf = luf_create_it();
    1.67 +      lpf->m = 0;
    1.68 +      lpf->B = NULL;
    1.69 +      lpf->n_max = 50;
    1.70 +      lpf->n = 0;
    1.71 +      lpf->R_ptr = lpf->R_len = NULL;
    1.72 +      lpf->S_ptr = lpf->S_len = NULL;
    1.73 +      lpf->scf = NULL;
    1.74 +      lpf->P_row = lpf->P_col = NULL;
    1.75 +      lpf->Q_row = lpf->Q_col = NULL;
    1.76 +      lpf->v_size = 1000;
    1.77 +      lpf->v_ptr = 0;
    1.78 +      lpf->v_ind = NULL;
    1.79 +      lpf->v_val = NULL;
    1.80 +      lpf->work1 = lpf->work2 = NULL;
    1.81 +      return lpf;
    1.82 +}
    1.83 +
    1.84 +/***********************************************************************
    1.85 +*  NAME
    1.86 +*
    1.87 +*  lpf_factorize - compute LP basis factorization
    1.88 +*
    1.89 +*  SYNOPSIS
    1.90 +*
    1.91 +*  #include "glplpf.h"
    1.92 +*  int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
    1.93 +*     (void *info, int j, int ind[], double val[]), void *info);
    1.94 +*
    1.95 +*  DESCRIPTION
    1.96 +*
    1.97 +*  The routine lpf_factorize computes the factorization of the basis
    1.98 +*  matrix B specified by the routine col.
    1.99 +*
   1.100 +*  The parameter lpf specified the basis factorization data structure
   1.101 +*  created with the routine lpf_create_it.
   1.102 +*
   1.103 +*  The parameter m specifies the order of B, m > 0.
   1.104 +*
   1.105 +*  The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
   1.106 +*  number of j-th column of B in some original matrix. The array bh is
   1.107 +*  optional and can be specified as NULL.
   1.108 +*
   1.109 +*  The formal routine col specifies the matrix B to be factorized. To
   1.110 +*  obtain j-th column of A the routine lpf_factorize calls the routine
   1.111 +*  col with the parameter j (1 <= j <= n). In response the routine col
   1.112 +*  should store row indices and numerical values of non-zero elements
   1.113 +*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
   1.114 +*  respectively, where len is the number of non-zeros in j-th column
   1.115 +*  returned on exit. Neither zero nor duplicate elements are allowed.
   1.116 +*
   1.117 +*  The parameter info is a transit pointer passed to the routine col.
   1.118 +*
   1.119 +*  RETURNS
   1.120 +*
   1.121 +*  0  The factorization has been successfully computed.
   1.122 +*
   1.123 +*  LPF_ESING
   1.124 +*     The specified matrix is singular within the working precision.
   1.125 +*
   1.126 +*  LPF_ECOND
   1.127 +*     The specified matrix is ill-conditioned.
   1.128 +*
   1.129 +*  For more details see comments to the routine luf_factorize. */
   1.130 +
   1.131 +int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
   1.132 +      (void *info, int j, int ind[], double val[]), void *info)
   1.133 +{     int k, ret;
   1.134 +#if _GLPLPF_DEBUG
   1.135 +      int i, j, len, *ind;
   1.136 +      double *B, *val;
   1.137 +#endif
   1.138 +      xassert(bh == bh);
   1.139 +      if (m < 1)
   1.140 +         xfault("lpf_factorize: m = %d; invalid parameter\n", m);
   1.141 +      if (m > M_MAX)
   1.142 +         xfault("lpf_factorize: m = %d; matrix too big\n", m);
   1.143 +      lpf->m0 = lpf->m = m;
   1.144 +      /* invalidate the factorization */
   1.145 +      lpf->valid = 0;
   1.146 +      /* allocate/reallocate arrays, if necessary */
   1.147 +      if (lpf->R_ptr == NULL)
   1.148 +         lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int));
   1.149 +      if (lpf->R_len == NULL)
   1.150 +         lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int));
   1.151 +      if (lpf->S_ptr == NULL)
   1.152 +         lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int));
   1.153 +      if (lpf->S_len == NULL)
   1.154 +         lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int));
   1.155 +      if (lpf->scf == NULL)
   1.156 +         lpf->scf = scf_create_it(lpf->n_max);
   1.157 +      if (lpf->v_ind == NULL)
   1.158 +         lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int));
   1.159 +      if (lpf->v_val == NULL)
   1.160 +         lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double));
   1.161 +      if (lpf->m0_max < m)
   1.162 +      {  if (lpf->P_row != NULL) xfree(lpf->P_row);
   1.163 +         if (lpf->P_col != NULL) xfree(lpf->P_col);
   1.164 +         if (lpf->Q_row != NULL) xfree(lpf->Q_row);
   1.165 +         if (lpf->Q_col != NULL) xfree(lpf->Q_col);
   1.166 +         if (lpf->work1 != NULL) xfree(lpf->work1);
   1.167 +         if (lpf->work2 != NULL) xfree(lpf->work2);
   1.168 +         lpf->m0_max = m + 100;
   1.169 +         lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
   1.170 +         lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
   1.171 +         lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
   1.172 +         lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
   1.173 +         lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
   1.174 +         lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
   1.175 +      }
   1.176 +      /* try to factorize the basis matrix */
   1.177 +      switch (luf_factorize(lpf->luf, m, col, info))
   1.178 +      {  case 0:
   1.179 +            break;
   1.180 +         case LUF_ESING:
   1.181 +            ret = LPF_ESING;
   1.182 +            goto done;
   1.183 +         case LUF_ECOND:
   1.184 +            ret = LPF_ECOND;
   1.185 +            goto done;
   1.186 +         default:
   1.187 +            xassert(lpf != lpf);
   1.188 +      }
   1.189 +      /* the basis matrix has been successfully factorized */
   1.190 +      lpf->valid = 1;
   1.191 +#if _GLPLPF_DEBUG
   1.192 +      /* store the basis matrix for debugging */
   1.193 +      if (lpf->B != NULL) xfree(lpf->B);
   1.194 +      xassert(m <= 32767);
   1.195 +      lpf->B = B = xcalloc(1+m*m, sizeof(double));
   1.196 +      ind = xcalloc(1+m, sizeof(int));
   1.197 +      val = xcalloc(1+m, sizeof(double));
   1.198 +      for (k = 1; k <= m * m; k++)
   1.199 +         B[k] = 0.0;
   1.200 +      for (j = 1; j <= m; j++)
   1.201 +      {  len = col(info, j, ind, val);
   1.202 +         xassert(0 <= len && len <= m);
   1.203 +         for (k = 1; k <= len; k++)
   1.204 +         {  i = ind[k];
   1.205 +            xassert(1 <= i && i <= m);
   1.206 +            xassert(B[(i - 1) * m + j] == 0.0);
   1.207 +            xassert(val[k] != 0.0);
   1.208 +            B[(i - 1) * m + j] = val[k];
   1.209 +         }
   1.210 +      }
   1.211 +      xfree(ind);
   1.212 +      xfree(val);
   1.213 +#endif
   1.214 +      /* B = B0, so there are no additional rows/columns */
   1.215 +      lpf->n = 0;
   1.216 +      /* reset the Schur complement factorization */
   1.217 +      scf_reset_it(lpf->scf);
   1.218 +      /* P := Q := I */
   1.219 +      for (k = 1; k <= m; k++)
   1.220 +      {  lpf->P_row[k] = lpf->P_col[k] = k;
   1.221 +         lpf->Q_row[k] = lpf->Q_col[k] = k;
   1.222 +      }
   1.223 +      /* make all SVA locations free */
   1.224 +      lpf->v_ptr = 1;
   1.225 +      ret = 0;
   1.226 +done: /* return to the calling program */
   1.227 +      return ret;
   1.228 +}
   1.229 +
   1.230 +/***********************************************************************
   1.231 +*  The routine r_prod computes the product y := y + alpha * R * x,
   1.232 +*  where x is a n-vector, alpha is a scalar, y is a m0-vector.
   1.233 +*
   1.234 +*  Since matrix R is available by columns, the product is computed as
   1.235 +*  a linear combination:
   1.236 +*
   1.237 +*     y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]),
   1.238 +*
   1.239 +*  where R[j] is j-th column of R. */
   1.240 +
   1.241 +static void r_prod(LPF *lpf, double y[], double a, const double x[])
   1.242 +{     int n = lpf->n;
   1.243 +      int *R_ptr = lpf->R_ptr;
   1.244 +      int *R_len = lpf->R_len;
   1.245 +      int *v_ind = lpf->v_ind;
   1.246 +      double *v_val = lpf->v_val;
   1.247 +      int j, beg, end, ptr;
   1.248 +      double t;
   1.249 +      for (j = 1; j <= n; j++)
   1.250 +      {  if (x[j] == 0.0) continue;
   1.251 +         /* y := y + alpha * R[j] * x[j] */
   1.252 +         t = a * x[j];
   1.253 +         beg = R_ptr[j];
   1.254 +         end = beg + R_len[j];
   1.255 +         for (ptr = beg; ptr < end; ptr++)
   1.256 +            y[v_ind[ptr]] += t * v_val[ptr];
   1.257 +      }
   1.258 +      return;
   1.259 +}
   1.260 +
   1.261 +/***********************************************************************
   1.262 +*  The routine rt_prod computes the product y := y + alpha * R' * x,
   1.263 +*  where R' is a matrix transposed to R, x is a m0-vector, alpha is a
   1.264 +*  scalar, y is a n-vector.
   1.265 +*
   1.266 +*  Since matrix R is available by columns, the product components are
   1.267 +*  computed as inner products:
   1.268 +*
   1.269 +*     y[j] := y[j] + alpha * (j-th column of R) * x
   1.270 +*
   1.271 +*  for j = 1, 2, ..., n. */
   1.272 +
   1.273 +static void rt_prod(LPF *lpf, double y[], double a, const double x[])
   1.274 +{     int n = lpf->n;
   1.275 +      int *R_ptr = lpf->R_ptr;
   1.276 +      int *R_len = lpf->R_len;
   1.277 +      int *v_ind = lpf->v_ind;
   1.278 +      double *v_val = lpf->v_val;
   1.279 +      int j, beg, end, ptr;
   1.280 +      double t;
   1.281 +      for (j = 1; j <= n; j++)
   1.282 +      {  /* t := (j-th column of R) * x */
   1.283 +         t = 0.0;
   1.284 +         beg = R_ptr[j];
   1.285 +         end = beg + R_len[j];
   1.286 +         for (ptr = beg; ptr < end; ptr++)
   1.287 +            t += v_val[ptr] * x[v_ind[ptr]];
   1.288 +         /* y[j] := y[j] + alpha * t */
   1.289 +         y[j] += a * t;
   1.290 +      }
   1.291 +      return;
   1.292 +}
   1.293 +
   1.294 +/***********************************************************************
   1.295 +*  The routine s_prod computes the product y := y + alpha * S * x,
   1.296 +*  where x is a m0-vector, alpha is a scalar, y is a n-vector.
   1.297 +*
   1.298 +*  Since matrix S is available by rows, the product components are
   1.299 +*  computed as inner products:
   1.300 +*
   1.301 +*     y[i] = y[i] + alpha * (i-th row of S) * x
   1.302 +*
   1.303 +*  for i = 1, 2, ..., n. */
   1.304 +
   1.305 +static void s_prod(LPF *lpf, double y[], double a, const double x[])
   1.306 +{     int n = lpf->n;
   1.307 +      int *S_ptr = lpf->S_ptr;
   1.308 +      int *S_len = lpf->S_len;
   1.309 +      int *v_ind = lpf->v_ind;
   1.310 +      double *v_val = lpf->v_val;
   1.311 +      int i, beg, end, ptr;
   1.312 +      double t;
   1.313 +      for (i = 1; i <= n; i++)
   1.314 +      {  /* t := (i-th row of S) * x */
   1.315 +         t = 0.0;
   1.316 +         beg = S_ptr[i];
   1.317 +         end = beg + S_len[i];
   1.318 +         for (ptr = beg; ptr < end; ptr++)
   1.319 +            t += v_val[ptr] * x[v_ind[ptr]];
   1.320 +         /* y[i] := y[i] + alpha * t */
   1.321 +         y[i] += a * t;
   1.322 +      }
   1.323 +      return;
   1.324 +}
   1.325 +
   1.326 +/***********************************************************************
   1.327 +*  The routine st_prod computes the product y := y + alpha * S' * x,
   1.328 +*  where S' is a matrix transposed to S, x is a n-vector, alpha is a
   1.329 +*  scalar, y is m0-vector.
   1.330 +*
   1.331 +*  Since matrix R is available by rows, the product is computed as a
   1.332 +*  linear combination:
   1.333 +*
   1.334 +*     y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]),
   1.335 +*
   1.336 +*  where S'[i] is i-th row of S. */
   1.337 +
   1.338 +static void st_prod(LPF *lpf, double y[], double a, const double x[])
   1.339 +{     int n = lpf->n;
   1.340 +      int *S_ptr = lpf->S_ptr;
   1.341 +      int *S_len = lpf->S_len;
   1.342 +      int *v_ind = lpf->v_ind;
   1.343 +      double *v_val = lpf->v_val;
   1.344 +      int i, beg, end, ptr;
   1.345 +      double t;
   1.346 +      for (i = 1; i <= n; i++)
   1.347 +      {  if (x[i] == 0.0) continue;
   1.348 +         /* y := y + alpha * S'[i] * x[i] */
   1.349 +         t = a * x[i];
   1.350 +         beg = S_ptr[i];
   1.351 +         end = beg + S_len[i];
   1.352 +         for (ptr = beg; ptr < end; ptr++)
   1.353 +            y[v_ind[ptr]] += t * v_val[ptr];
   1.354 +      }
   1.355 +      return;
   1.356 +}
   1.357 +
   1.358 +#if _GLPLPF_DEBUG
   1.359 +/***********************************************************************
   1.360 +*  The routine check_error computes the maximal relative error between
   1.361 +*  left- and right-hand sides for the system B * x = b (if tr is zero)
   1.362 +*  or B' * x = b (if tr is non-zero), where B' is a matrix transposed
   1.363 +*  to B. (This routine is intended for debugging only.) */
   1.364 +
   1.365 +static void check_error(LPF *lpf, int tr, const double x[],
   1.366 +      const double b[])
   1.367 +{     int m = lpf->m;
   1.368 +      double *B = lpf->B;
   1.369 +      int i, j;
   1.370 +      double  d, dmax = 0.0, s, t, tmax;
   1.371 +      for (i = 1; i <= m; i++)
   1.372 +      {  s = 0.0;
   1.373 +         tmax = 1.0;
   1.374 +         for (j = 1; j <= m; j++)
   1.375 +         {  if (!tr)
   1.376 +               t = B[m * (i - 1) + j] * x[j];
   1.377 +            else
   1.378 +               t = B[m * (j - 1) + i] * x[j];
   1.379 +            if (tmax < fabs(t)) tmax = fabs(t);
   1.380 +            s += t;
   1.381 +         }
   1.382 +         d = fabs(s - b[i]) / tmax;
   1.383 +         if (dmax < d) dmax = d;
   1.384 +      }
   1.385 +      if (dmax > 1e-8)
   1.386 +         xprintf("%s: dmax = %g; relative error too large\n",
   1.387 +            !tr ? "lpf_ftran" : "lpf_btran", dmax);
   1.388 +      return;
   1.389 +}
   1.390 +#endif
   1.391 +
   1.392 +/***********************************************************************
   1.393 +*  NAME
   1.394 +*
   1.395 +*  lpf_ftran - perform forward transformation (solve system B*x = b)
   1.396 +*
   1.397 +*  SYNOPSIS
   1.398 +*
   1.399 +*  #include "glplpf.h"
   1.400 +*  void lpf_ftran(LPF *lpf, double x[]);
   1.401 +*
   1.402 +*  DESCRIPTION
   1.403 +*
   1.404 +*  The routine lpf_ftran performs forward transformation, i.e. solves
   1.405 +*  the system B*x = b, where B is the basis matrix, x is the vector of
   1.406 +*  unknowns to be computed, b is the vector of right-hand sides.
   1.407 +*
   1.408 +*  On entry elements of the vector b should be stored in dense format
   1.409 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.410 +*  the routine stores elements of the vector x in the same locations.
   1.411 +*
   1.412 +*  BACKGROUND
   1.413 +*
   1.414 +*  Solution of the system B * x = b can be obtained by solving the
   1.415 +*  following augmented system:
   1.416 +*
   1.417 +*     ( B  F^) ( x )   ( b )
   1.418 +*     (      ) (   ) = (   )
   1.419 +*     ( G^ H^) ( y )   ( 0 )
   1.420 +*
   1.421 +*  which, using the main equality, can be written as follows:
   1.422 +*
   1.423 +*       ( L0 0 ) ( U0 R )   ( x )   ( b )
   1.424 +*     P (      ) (      ) Q (   ) = (   )
   1.425 +*       ( S  I ) ( 0  C )   ( y )   ( 0 )
   1.426 +*
   1.427 +*  therefore,
   1.428 +*
   1.429 +*     ( x )      ( U0 R )-1 ( L0 0 )-1    ( b )
   1.430 +*     (   ) = Q' (      )   (      )   P' (   )
   1.431 +*     ( y )      ( 0  C )   ( S  I )      ( 0 )
   1.432 +*
   1.433 +*  Thus, computing the solution includes the following steps:
   1.434 +*
   1.435 +*  1. Compute
   1.436 +*
   1.437 +*     ( f )      ( b )
   1.438 +*     (   ) = P' (   )
   1.439 +*     ( g )      ( 0 )
   1.440 +*
   1.441 +*  2. Solve the system
   1.442 +*
   1.443 +*     ( f1 )   ( L0 0 )-1 ( f )      ( L0 0 ) ( f1 )   ( f )
   1.444 +*     (    ) = (      )   (   )  =>  (      ) (    ) = (   )
   1.445 +*     ( g1 )   ( S  I )   ( g )      ( S  I ) ( g1 )   ( g )
   1.446 +*
   1.447 +*     from which it follows that:
   1.448 +*
   1.449 +*     { L0 * f1      = f      f1 = inv(L0) * f
   1.450 +*     {                   =>
   1.451 +*     {  S * f1 + g1 = g      g1 = g - S * f1
   1.452 +*
   1.453 +*  3. Solve the system
   1.454 +*
   1.455 +*     ( f2 )   ( U0 R )-1 ( f1 )      ( U0 R ) ( f2 )   ( f1 )
   1.456 +*     (    ) = (      )   (    )  =>  (      ) (    ) = (    )
   1.457 +*     ( g2 )   ( 0  C )   ( g1 )      ( 0  C ) ( g2 )   ( g1 )
   1.458 +*
   1.459 +*     from which it follows that:
   1.460 +*
   1.461 +*     { U0 * f2 + R * g2 = f1      f2 = inv(U0) * (f1 - R * g2)
   1.462 +*     {                        =>
   1.463 +*     {           C * g2 = g1      g2 = inv(C) * g1
   1.464 +*
   1.465 +*  4. Compute
   1.466 +*
   1.467 +*     ( x )      ( f2 )
   1.468 +*     (   ) = Q' (    )
   1.469 +*     ( y )      ( g2 )                                               */
   1.470 +
   1.471 +void lpf_ftran(LPF *lpf, double x[])
   1.472 +{     int m0 = lpf->m0;
   1.473 +      int m = lpf->m;
   1.474 +      int n  = lpf->n;
   1.475 +      int *P_col = lpf->P_col;
   1.476 +      int *Q_col = lpf->Q_col;
   1.477 +      double *fg = lpf->work1;
   1.478 +      double *f = fg;
   1.479 +      double *g = fg + m0;
   1.480 +      int i, ii;
   1.481 +#if _GLPLPF_DEBUG
   1.482 +      double *b;
   1.483 +#endif
   1.484 +      if (!lpf->valid)
   1.485 +         xfault("lpf_ftran: the factorization is not valid\n");
   1.486 +      xassert(0 <= m && m <= m0 + n);
   1.487 +#if _GLPLPF_DEBUG
   1.488 +      /* save the right-hand side vector */
   1.489 +      b = xcalloc(1+m, sizeof(double));
   1.490 +      for (i = 1; i <= m; i++) b[i] = x[i];
   1.491 +#endif
   1.492 +      /* (f g) := inv(P) * (b 0) */
   1.493 +      for (i = 1; i <= m0 + n; i++)
   1.494 +         fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0);
   1.495 +      /* f1 := inv(L0) * f */
   1.496 +      luf_f_solve(lpf->luf, 0, f);
   1.497 +      /* g1 := g - S * f1 */
   1.498 +      s_prod(lpf, g, -1.0, f);
   1.499 +      /* g2 := inv(C) * g1 */
   1.500 +      scf_solve_it(lpf->scf, 0, g);
   1.501 +      /* f2 := inv(U0) * (f1 - R * g2) */
   1.502 +      r_prod(lpf, f, -1.0, g);
   1.503 +      luf_v_solve(lpf->luf, 0, f);
   1.504 +      /* (x y) := inv(Q) * (f2 g2) */
   1.505 +      for (i = 1; i <= m; i++)
   1.506 +         x[i] = fg[Q_col[i]];
   1.507 +#if _GLPLPF_DEBUG
   1.508 +      /* check relative error in solution */
   1.509 +      check_error(lpf, 0, x, b);
   1.510 +      xfree(b);
   1.511 +#endif
   1.512 +      return;
   1.513 +}
   1.514 +
   1.515 +/***********************************************************************
   1.516 +*  NAME
   1.517 +*
   1.518 +*  lpf_btran - perform backward transformation (solve system B'*x = b)
   1.519 +*
   1.520 +*  SYNOPSIS
   1.521 +*
   1.522 +*  #include "glplpf.h"
   1.523 +*  void lpf_btran(LPF *lpf, double x[]);
   1.524 +*
   1.525 +*  DESCRIPTION
   1.526 +*
   1.527 +*  The routine lpf_btran performs backward transformation, i.e. solves
   1.528 +*  the system B'*x = b, where B' is a matrix transposed to the basis
   1.529 +*  matrix B, x is the vector of unknowns to be computed, b is the vector
   1.530 +*  of right-hand sides.
   1.531 +*
   1.532 +*  On entry elements of the vector b should be stored in dense format
   1.533 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.534 +*  the routine stores elements of the vector x in the same locations.
   1.535 +*
   1.536 +*  BACKGROUND
   1.537 +*
   1.538 +*  Solution of the system B' * x = b, where B' is a matrix transposed
   1.539 +*  to B, can be obtained by solving the following augmented system:
   1.540 +*
   1.541 +*     ( B  F^)T ( x )   ( b )
   1.542 +*     (      )  (   ) = (   )
   1.543 +*     ( G^ H^)  ( y )   ( 0 )
   1.544 +*
   1.545 +*  which, using the main equality, can be written as follows:
   1.546 +*
   1.547 +*      T ( U0 R )T ( L0 0 )T  T ( x )   ( b )
   1.548 +*     Q  (      )  (      )  P  (   ) = (   )
   1.549 +*        ( 0  C )  ( S  I )     ( y )   ( 0 )
   1.550 +*
   1.551 +*  or, equivalently, as follows:
   1.552 +*
   1.553 +*        ( U'0 0 ) ( L'0 S')    ( x )   ( b )
   1.554 +*     Q' (       ) (       ) P' (   ) = (   )
   1.555 +*        ( R'  C') ( 0   I )    ( y )   ( 0 )
   1.556 +*
   1.557 +*  therefore,
   1.558 +*
   1.559 +*     ( x )     ( L'0 S')-1 ( U'0 0 )-1   ( b )
   1.560 +*     (   ) = P (       )   (       )   Q (   )
   1.561 +*     ( y )     ( 0   I )   ( R'  C')     ( 0 )
   1.562 +*
   1.563 +*  Thus, computing the solution includes the following steps:
   1.564 +*
   1.565 +*  1. Compute
   1.566 +*
   1.567 +*     ( f )     ( b )
   1.568 +*     (   ) = Q (   )
   1.569 +*     ( g )     ( 0 )
   1.570 +*
   1.571 +*  2. Solve the system
   1.572 +*
   1.573 +*     ( f1 )   ( U'0 0 )-1 ( f )      ( U'0 0 ) ( f1 )   ( f )
   1.574 +*     (    ) = (       )   (   )  =>  (       ) (    ) = (   )
   1.575 +*     ( g1 )   ( R'  C')   ( g )      ( R'  C') ( g1 )   ( g )
   1.576 +*
   1.577 +*     from which it follows that:
   1.578 +*
   1.579 +*     { U'0 * f1           = f      f1 = inv(U'0) * f
   1.580 +*     {                         =>
   1.581 +*     { R'  * f1 + C' * g1 = g      g1 = inv(C') * (g - R' * f1)
   1.582 +*
   1.583 +*  3. Solve the system
   1.584 +*
   1.585 +*     ( f2 )   ( L'0 S')-1 ( f1 )      ( L'0 S') ( f2 )   ( f1 )
   1.586 +*     (    ) = (       )   (    )  =>  (       ) (    ) = (    )
   1.587 +*     ( g2 )   ( 0   I )   ( g1 )      ( 0   I ) ( g2 )   ( g1 )
   1.588 +*
   1.589 +*     from which it follows that:
   1.590 +*
   1.591 +*     { L'0 * f2 + S' * g2 = f1
   1.592 +*     {                          =>  f2 = inv(L'0) * ( f1 - S' * g2)
   1.593 +*     {                 g2 = g1
   1.594 +*
   1.595 +*  4. Compute
   1.596 +*
   1.597 +*     ( x )     ( f2 )
   1.598 +*     (   ) = P (    )
   1.599 +*     ( y )     ( g2 )                                                */
   1.600 +
   1.601 +void lpf_btran(LPF *lpf, double x[])
   1.602 +{     int m0 = lpf->m0;
   1.603 +      int m = lpf->m;
   1.604 +      int n = lpf->n;
   1.605 +      int *P_row = lpf->P_row;
   1.606 +      int *Q_row = lpf->Q_row;
   1.607 +      double *fg = lpf->work1;
   1.608 +      double *f = fg;
   1.609 +      double *g = fg + m0;
   1.610 +      int i, ii;
   1.611 +#if _GLPLPF_DEBUG
   1.612 +      double *b;
   1.613 +#endif
   1.614 +      if (!lpf->valid)
   1.615 +         xfault("lpf_btran: the factorization is not valid\n");
   1.616 +      xassert(0 <= m && m <= m0 + n);
   1.617 +#if _GLPLPF_DEBUG
   1.618 +      /* save the right-hand side vector */
   1.619 +      b = xcalloc(1+m, sizeof(double));
   1.620 +      for (i = 1; i <= m; i++) b[i] = x[i];
   1.621 +#endif
   1.622 +      /* (f g) := Q * (b 0) */
   1.623 +      for (i = 1; i <= m0 + n; i++)
   1.624 +         fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0);
   1.625 +      /* f1 := inv(U'0) * f */
   1.626 +      luf_v_solve(lpf->luf, 1, f);
   1.627 +      /* g1 := inv(C') * (g - R' * f1) */
   1.628 +      rt_prod(lpf, g, -1.0, f);
   1.629 +      scf_solve_it(lpf->scf, 1, g);
   1.630 +      /* g2 := g1 */
   1.631 +      g = g;
   1.632 +      /* f2 := inv(L'0) * (f1 - S' * g2) */
   1.633 +      st_prod(lpf, f, -1.0, g);
   1.634 +      luf_f_solve(lpf->luf, 1, f);
   1.635 +      /* (x y) := P * (f2 g2) */
   1.636 +      for (i = 1; i <= m; i++)
   1.637 +         x[i] = fg[P_row[i]];
   1.638 +#if _GLPLPF_DEBUG
   1.639 +      /* check relative error in solution */
   1.640 +      check_error(lpf, 1, x, b);
   1.641 +      xfree(b);
   1.642 +#endif
   1.643 +      return;
   1.644 +}
   1.645 +
   1.646 +/***********************************************************************
   1.647 +*  The routine enlarge_sva enlarges the Sparse Vector Area to new_size
   1.648 +*  locations by reallocating the arrays v_ind and v_val. */
   1.649 +
   1.650 +static void enlarge_sva(LPF *lpf, int new_size)
   1.651 +{     int v_size = lpf->v_size;
   1.652 +      int used = lpf->v_ptr - 1;
   1.653 +      int *v_ind = lpf->v_ind;
   1.654 +      double *v_val = lpf->v_val;
   1.655 +      xassert(v_size < new_size);
   1.656 +      while (v_size < new_size) v_size += v_size;
   1.657 +      lpf->v_size = v_size;
   1.658 +      lpf->v_ind = xcalloc(1+v_size, sizeof(int));
   1.659 +      lpf->v_val = xcalloc(1+v_size, sizeof(double));
   1.660 +      xassert(used >= 0);
   1.661 +      memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int));
   1.662 +      memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double));
   1.663 +      xfree(v_ind);
   1.664 +      xfree(v_val);
   1.665 +      return;
   1.666 +}
   1.667 +
   1.668 +/***********************************************************************
   1.669 +*  NAME
   1.670 +*
   1.671 +*  lpf_update_it - update LP basis factorization
   1.672 +*
   1.673 +*  SYNOPSIS
   1.674 +*
   1.675 +*  #include "glplpf.h"
   1.676 +*  int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
   1.677 +*     const double val[]);
   1.678 +*
   1.679 +*  DESCRIPTION
   1.680 +*
   1.681 +*  The routine lpf_update_it updates the factorization of the basis
   1.682 +*  matrix B after replacing its j-th column by a new vector.
   1.683 +*
   1.684 +*  The parameter j specifies the number of column of B, which has been
   1.685 +*  replaced, 1 <= j <= m, where m is the order of B.
   1.686 +*
   1.687 +*  The parameter bh specifies the basis header entry for the new column
   1.688 +*  of B, which is the number of the new column in some original matrix.
   1.689 +*  This parameter is optional and can be specified as 0.
   1.690 +*
   1.691 +*  Row indices and numerical values of non-zero elements of the new
   1.692 +*  column of B should be placed in locations ind[1], ..., ind[len] and
   1.693 +*  val[1], ..., val[len], resp., where len is the number of non-zeros
   1.694 +*  in the column. Neither zero nor duplicate elements are allowed.
   1.695 +*
   1.696 +*  RETURNS
   1.697 +*
   1.698 +*  0  The factorization has been successfully updated.
   1.699 +*
   1.700 +*  LPF_ESING
   1.701 +*     New basis B is singular within the working precision.
   1.702 +*
   1.703 +*  LPF_ELIMIT
   1.704 +*     Maximal number of additional rows and columns has been reached.
   1.705 +*
   1.706 +*  BACKGROUND
   1.707 +*
   1.708 +*  Let j-th column of the current basis matrix B have to be replaced by
   1.709 +*  a new column a. This replacement is equivalent to removing the old
   1.710 +*  j-th column by fixing it at zero and introducing the new column as
   1.711 +*  follows:
   1.712 +*
   1.713 +*                   ( B   F^| a )
   1.714 +*     ( B  F^)      (       |   )
   1.715 +*     (      ) ---> ( G^  H^| 0 )
   1.716 +*     ( G^ H^)      (-------+---)
   1.717 +*                   ( e'j 0 | 0 )
   1.718 +*
   1.719 +*  where ej is a unit vector with 1 in j-th position which used to fix
   1.720 +*  the old j-th column of B (at zero). Then using the main equality we
   1.721 +*  have:
   1.722 +*
   1.723 +*     ( B   F^| a )            ( B0  F | f )
   1.724 +*     (       |   )   ( P  0 ) (       |   ) ( Q  0 )
   1.725 +*     ( G^  H^| 0 ) = (      ) ( G   H | g ) (      ) =
   1.726 +*     (-------+---)   ( 0  1 ) (-------+---) ( 0  1 )
   1.727 +*     ( e'j 0 | 0 )            ( v'  w'| 0 )
   1.728 +*
   1.729 +*       [   ( B0  F )|   ( f ) ]            [   ( B0 F )  |   ( f ) ]
   1.730 +*       [ P (       )| P (   ) ] ( Q  0 )   [ P (      ) Q| P (   ) ]
   1.731 +*     = [   ( G   H )|   ( g ) ] (      ) = [   ( G  H )  |   ( g ) ]
   1.732 +*       [------------+-------- ] ( 0  1 )   [-------------+---------]
   1.733 +*       [   ( v'  w')|     0   ]            [   ( v' w') Q|    0    ]
   1.734 +*
   1.735 +*  where:
   1.736 +*
   1.737 +*     ( a )     ( f )      ( f )        ( a )
   1.738 +*     (   ) = P (   )  =>  (   ) = P' * (   )
   1.739 +*     ( 0 )     ( g )      ( g )        ( 0 )
   1.740 +*
   1.741 +*                                 ( ej )      ( v )    ( v )     ( ej )
   1.742 +*     ( e'j  0 ) = ( v' w' ) Q => (    ) = Q' (   ) => (   ) = Q (    )
   1.743 +*                                 ( 0  )      ( w )    ( w )     ( 0  )
   1.744 +*
   1.745 +*  On the other hand:
   1.746 +*
   1.747 +*              ( B0| F  f )
   1.748 +*     ( P  0 ) (---+------) ( Q  0 )         ( B0    new F )
   1.749 +*     (      ) ( G | H  g ) (      ) = new P (             ) new Q
   1.750 +*     ( 0  1 ) (   |      ) ( 0  1 )         ( new G new H )
   1.751 +*              ( v'| w' 0 )
   1.752 +*
   1.753 +*  where:
   1.754 +*                               ( G )           ( H  g )
   1.755 +*     new F = ( F  f ), new G = (   ),  new H = (      ),
   1.756 +*                               ( v')           ( w' 0 )
   1.757 +*
   1.758 +*             ( P  0 )           ( Q  0 )
   1.759 +*     new P = (      ) , new Q = (      ) .
   1.760 +*             ( 0  1 )           ( 0  1 )
   1.761 +*
   1.762 +*  The factorization structure for the new augmented matrix remains the
   1.763 +*  same, therefore:
   1.764 +*
   1.765 +*           ( B0    new F )         ( L0     0 ) ( U0 new R )
   1.766 +*     new P (             ) new Q = (          ) (          )
   1.767 +*           ( new G new H )         ( new S  I ) ( 0  new C )
   1.768 +*
   1.769 +*  where:
   1.770 +*
   1.771 +*     new F = L0 * new R  =>
   1.772 +*
   1.773 +*     new R = inv(L0) * new F = inv(L0) * (F  f) = ( R  inv(L0)*f )
   1.774 +*
   1.775 +*     new G = new S * U0  =>
   1.776 +*
   1.777 +*                               ( G )             (     S      )
   1.778 +*     new S = new G * inv(U0) = (   ) * inv(U0) = (            )
   1.779 +*                               ( v')             ( v'*inv(U0) )
   1.780 +*
   1.781 +*     new H = new S * new R + new C  =>
   1.782 +*
   1.783 +*     new C = new H - new S * new R =
   1.784 +*
   1.785 +*             ( H  g )   (     S      )
   1.786 +*           = (      ) - (            ) * ( R  inv(L0)*f ) =
   1.787 +*             ( w' 0 )   ( v'*inv(U0) )
   1.788 +*
   1.789 +*             ( H - S*R           g - S*inv(L0)*f      )   ( C  x )
   1.790 +*           = (                                        ) = (      )
   1.791 +*             ( w'- v'*inv(U0)*R  -v'*inv(U0)*inv(L0)*f)   ( y' z )
   1.792 +*
   1.793 +*  Note that new C is resulted by expanding old C with new column x,
   1.794 +*  row y', and diagonal element z, where:
   1.795 +*
   1.796 +*     x = g - S * inv(L0) * f = g - S * (new column of R)
   1.797 +*
   1.798 +*     y = w - R'* inv(U'0)* v = w - R'* (new row of S)
   1.799 +*
   1.800 +*     z = - (new row of S) * (new column of R)
   1.801 +*
   1.802 +*  Finally, to replace old B by new B we have to permute j-th and last
   1.803 +*  (just added) columns of the matrix
   1.804 +*
   1.805 +*     ( B   F^| a )
   1.806 +*     (       |   )
   1.807 +*     ( G^  H^| 0 )
   1.808 +*     (-------+---)
   1.809 +*     ( e'j 0 | 0 )
   1.810 +*
   1.811 +*  and to keep the main equality do the same for matrix Q. */
   1.812 +
   1.813 +int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
   1.814 +      const double val[])
   1.815 +{     int m0 = lpf->m0;
   1.816 +      int m = lpf->m;
   1.817 +#if _GLPLPF_DEBUG
   1.818 +      double *B = lpf->B;
   1.819 +#endif
   1.820 +      int n = lpf->n;
   1.821 +      int *R_ptr = lpf->R_ptr;
   1.822 +      int *R_len = lpf->R_len;
   1.823 +      int *S_ptr = lpf->S_ptr;
   1.824 +      int *S_len = lpf->S_len;
   1.825 +      int *P_row = lpf->P_row;
   1.826 +      int *P_col = lpf->P_col;
   1.827 +      int *Q_row = lpf->Q_row;
   1.828 +      int *Q_col = lpf->Q_col;
   1.829 +      int v_ptr = lpf->v_ptr;
   1.830 +      int *v_ind = lpf->v_ind;
   1.831 +      double *v_val = lpf->v_val;
   1.832 +      double *a = lpf->work2; /* new column */
   1.833 +      double *fg = lpf->work1, *f = fg, *g = fg + m0;
   1.834 +      double *vw = lpf->work2, *v = vw, *w = vw + m0;
   1.835 +      double *x = g, *y = w, z;
   1.836 +      int i, ii, k, ret;
   1.837 +      xassert(bh == bh);
   1.838 +      if (!lpf->valid)
   1.839 +         xfault("lpf_update_it: the factorization is not valid\n");
   1.840 +      if (!(1 <= j && j <= m))
   1.841 +         xfault("lpf_update_it: j = %d; column number out of range\n",
   1.842 +            j);
   1.843 +      xassert(0 <= m && m <= m0 + n);
   1.844 +      /* check if the basis factorization can be expanded */
   1.845 +      if (n == lpf->n_max)
   1.846 +      {  lpf->valid = 0;
   1.847 +         ret = LPF_ELIMIT;
   1.848 +         goto done;
   1.849 +      }
   1.850 +      /* convert new j-th column of B to dense format */
   1.851 +      for (i = 1; i <= m; i++)
   1.852 +         a[i] = 0.0;
   1.853 +      for (k = 1; k <= len; k++)
   1.854 +      {  i = ind[k];
   1.855 +         if (!(1 <= i && i <= m))
   1.856 +            xfault("lpf_update_it: ind[%d] = %d; row number out of rang"
   1.857 +               "e\n", k, i);
   1.858 +         if (a[i] != 0.0)
   1.859 +            xfault("lpf_update_it: ind[%d] = %d; duplicate row index no"
   1.860 +               "t allowed\n", k, i);
   1.861 +         if (val[k] == 0.0)
   1.862 +            xfault("lpf_update_it: val[%d] = %g; zero element not allow"
   1.863 +               "ed\n", k, val[k]);
   1.864 +         a[i] = val[k];
   1.865 +      }
   1.866 +#if _GLPLPF_DEBUG
   1.867 +      /* change column in the basis matrix for debugging */
   1.868 +      for (i = 1; i <= m; i++)
   1.869 +         B[(i - 1) * m + j] = a[i];
   1.870 +#endif
   1.871 +      /* (f g) := inv(P) * (a 0) */
   1.872 +      for (i = 1; i <= m0+n; i++)
   1.873 +         fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0);
   1.874 +      /* (v w) := Q * (ej 0) */
   1.875 +      for (i = 1; i <= m0+n; i++) vw[i] = 0.0;
   1.876 +      vw[Q_col[j]] = 1.0;
   1.877 +      /* f1 := inv(L0) * f (new column of R) */
   1.878 +      luf_f_solve(lpf->luf, 0, f);
   1.879 +      /* v1 := inv(U'0) * v (new row of S) */
   1.880 +      luf_v_solve(lpf->luf, 1, v);
   1.881 +      /* we need at most 2 * m0 available locations in the SVA to store
   1.882 +         new column of matrix R and new row of matrix S */
   1.883 +      if (lpf->v_size < v_ptr + m0 + m0)
   1.884 +      {  enlarge_sva(lpf, v_ptr + m0 + m0);
   1.885 +         v_ind = lpf->v_ind;
   1.886 +         v_val = lpf->v_val;
   1.887 +      }
   1.888 +      /* store new column of R */
   1.889 +      R_ptr[n+1] = v_ptr;
   1.890 +      for (i = 1; i <= m0; i++)
   1.891 +      {  if (f[i] != 0.0)
   1.892 +            v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++;
   1.893 +      }
   1.894 +      R_len[n+1] = v_ptr - lpf->v_ptr;
   1.895 +      lpf->v_ptr = v_ptr;
   1.896 +      /* store new row of S */
   1.897 +      S_ptr[n+1] = v_ptr;
   1.898 +      for (i = 1; i <= m0; i++)
   1.899 +      {  if (v[i] != 0.0)
   1.900 +            v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++;
   1.901 +      }
   1.902 +      S_len[n+1] = v_ptr - lpf->v_ptr;
   1.903 +      lpf->v_ptr = v_ptr;
   1.904 +      /* x := g - S * f1 (new column of C) */
   1.905 +      s_prod(lpf, x, -1.0, f);
   1.906 +      /* y := w - R' * v1 (new row of C) */
   1.907 +      rt_prod(lpf, y, -1.0, v);
   1.908 +      /* z := - v1 * f1 (new diagonal element of C) */
   1.909 +      z = 0.0;
   1.910 +      for (i = 1; i <= m0; i++) z -= v[i] * f[i];
   1.911 +      /* update factorization of new matrix C */
   1.912 +      switch (scf_update_exp(lpf->scf, x, y, z))
   1.913 +      {  case 0:
   1.914 +            break;
   1.915 +         case SCF_ESING:
   1.916 +            lpf->valid = 0;
   1.917 +            ret = LPF_ESING;
   1.918 +            goto done;
   1.919 +         case SCF_ELIMIT:
   1.920 +            xassert(lpf != lpf);
   1.921 +         default:
   1.922 +            xassert(lpf != lpf);
   1.923 +      }
   1.924 +      /* expand matrix P */
   1.925 +      P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1;
   1.926 +      /* expand matrix Q */
   1.927 +      Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1;
   1.928 +      /* permute j-th and last (just added) column of matrix Q */
   1.929 +      i = Q_col[j], ii = Q_col[m0+n+1];
   1.930 +      Q_row[i] = m0+n+1, Q_col[m0+n+1] = i;
   1.931 +      Q_row[ii] = j, Q_col[j] = ii;
   1.932 +      /* increase the number of additional rows and columns */
   1.933 +      lpf->n++;
   1.934 +      xassert(lpf->n <= lpf->n_max);
   1.935 +      /* the factorization has been successfully updated */
   1.936 +      ret = 0;
   1.937 +done: /* return to the calling program */
   1.938 +      return ret;
   1.939 +}
   1.940 +
   1.941 +/***********************************************************************
   1.942 +*  NAME
   1.943 +*
   1.944 +*  lpf_delete_it - delete LP basis factorization
   1.945 +*
   1.946 +*  SYNOPSIS
   1.947 +*
   1.948 +*  #include "glplpf.h"
   1.949 +*  void lpf_delete_it(LPF *lpf)
   1.950 +*
   1.951 +*  DESCRIPTION
   1.952 +*
   1.953 +*  The routine lpf_delete_it deletes LP basis factorization specified
   1.954 +*  by the parameter lpf and frees all memory allocated to this program
   1.955 +*  object. */
   1.956 +
   1.957 +void lpf_delete_it(LPF *lpf)
   1.958 +{     luf_delete_it(lpf->luf);
   1.959 +#if _GLPLPF_DEBUG
   1.960 +      if (lpf->B != NULL) xfree(lpf->B);
   1.961 +#else
   1.962 +      xassert(lpf->B == NULL);
   1.963 +#endif
   1.964 +      if (lpf->R_ptr != NULL) xfree(lpf->R_ptr);
   1.965 +      if (lpf->R_len != NULL) xfree(lpf->R_len);
   1.966 +      if (lpf->S_ptr != NULL) xfree(lpf->S_ptr);
   1.967 +      if (lpf->S_len != NULL) xfree(lpf->S_len);
   1.968 +      if (lpf->scf != NULL) scf_delete_it(lpf->scf);
   1.969 +      if (lpf->P_row != NULL) xfree(lpf->P_row);
   1.970 +      if (lpf->P_col != NULL) xfree(lpf->P_col);
   1.971 +      if (lpf->Q_row != NULL) xfree(lpf->Q_row);
   1.972 +      if (lpf->Q_col != NULL) xfree(lpf->Q_col);
   1.973 +      if (lpf->v_ind != NULL) xfree(lpf->v_ind);
   1.974 +      if (lpf->v_val != NULL) xfree(lpf->v_val);
   1.975 +      if (lpf->work1 != NULL) xfree(lpf->work1);
   1.976 +      if (lpf->work2 != NULL) xfree(lpf->work2);
   1.977 +      xfree(lpf);
   1.978 +      return;
   1.979 +}
   1.980 +
   1.981 +/* eof */