praat/external/glpk/glplpf.c

/* glplpf.c (LP basis factorization, Schur complement version) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
*  E-mail: <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glplpf.h"
#include "glpenv.h"
#define xfault xerror

#define _GLPLPF_DEBUG 0

/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */

#define M_MAX 100000000 /* = 100*10^6 */
/* maximal order of the basis matrix */

/***********************************************************************
*  NAME
*
*  lpf_create_it - create LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  LPF *lpf_create_it(void);
*
*  DESCRIPTION
*
*  The routine lpf_create_it creates a program object, which represents
*  a factorization of LP basis.
*
*  RETURNS
*
*  The routine lpf_create_it returns a pointer to the object created. */

LPF *lpf_create_it(void)
{     LPF *lpf;
#if _GLPLPF_DEBUG
      xprintf("lpf_create_it: warning: debug mode enabled\n");
#endif
      lpf = xmalloc(sizeof(LPF));
      lpf->valid = 0;
      lpf->m0_max = lpf->m0 = 0;
      lpf->luf = luf_create_it();
      lpf->m = 0;
      lpf->B = NULL;
      lpf->n_max = 50;
      lpf->n = 0;
      lpf->R_ptr = lpf->R_len = NULL;
      lpf->S_ptr = lpf->S_len = NULL;
      lpf->scf = NULL;
      lpf->P_row = lpf->P_col = NULL;
      lpf->Q_row = lpf->Q_col = NULL;
      lpf->v_size = 1000;
      lpf->v_ptr = 0;
      lpf->v_ind = NULL;
      lpf->v_val = NULL;
      lpf->work1 = lpf->work2 = NULL;
      return lpf;
}

/***********************************************************************
*  NAME
*
*  lpf_factorize - compute LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
*     (void *info, int j, int ind[], double val[]), void *info);
*
*  DESCRIPTION
*
*  The routine lpf_factorize computes the factorization of the basis
*  matrix B specified by the routine col.
*
*  The parameter lpf specified the basis factorization data structure
*  created with the routine lpf_create_it.
*
*  The parameter m specifies the order of B, m > 0.
*
*  The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
*  number of j-th column of B in some original matrix. The array bh is
*  optional and can be specified as NULL.
*
*  The formal routine col specifies the matrix B to be factorized. To
*  obtain j-th column of A the routine lpf_factorize calls the routine
*  col with the parameter j (1 <= j <= n). In response the routine col
*  should store row indices and numerical values of non-zero elements
*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
*  respectively, where len is the number of non-zeros in j-th column
*  returned on exit. Neither zero nor duplicate elements are allowed.
*
*  The parameter info is a transit pointer passed to the routine col.
*
*  RETURNS
*
*  0  The factorization has been successfully computed.
*
*  LPF_ESING
*     The specified matrix is singular within the working precision.
*
*  LPF_ECOND
*     The specified matrix is ill-conditioned.
*
*  For more details see comments to the routine luf_factorize. */

int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
      (void *info, int j, int ind[], double val[]), void *info)
{     int k, ret;
#if _GLPLPF_DEBUG
      int i, j, len, *ind;
      double *B, *val;
#endif
      xassert(bh == bh);
      if (m < 1)
         xfault("lpf_factorize: m = %d; invalid parameter\n", m);
      if (m > M_MAX)
         xfault("lpf_factorize: m = %d; matrix too big\n", m);
      lpf->m0 = lpf->m = m;
      /* invalidate the factorization */
      lpf->valid = 0;
      /* allocate/reallocate arrays, if necessary */
      if (lpf->R_ptr == NULL)
         lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->R_len == NULL)
         lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->S_ptr == NULL)
         lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->S_len == NULL)
         lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->scf == NULL)
         lpf->scf = scf_create_it(lpf->n_max);
      if (lpf->v_ind == NULL)
         lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int));
      if (lpf->v_val == NULL)
         lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double));
      if (lpf->m0_max < m)
      {  if (lpf->P_row != NULL) xfree(lpf->P_row);
         if (lpf->P_col != NULL) xfree(lpf->P_col);
         if (lpf->Q_row != NULL) xfree(lpf->Q_row);
         if (lpf->Q_col != NULL) xfree(lpf->Q_col);
         if (lpf->work1 != NULL) xfree(lpf->work1);
         if (lpf->work2 != NULL) xfree(lpf->work2);
         lpf->m0_max = m + 100;
         lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
         lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
      }
      /* try to factorize the basis matrix */
      switch (luf_factorize(lpf->luf, m, col, info))
      {  case 0:
            break;
         case LUF_ESING:
            ret = LPF_ESING;
            goto done;
         case LUF_ECOND:
            ret = LPF_ECOND;
            goto done;
         default:
            xassert(lpf != lpf);
      }
      /* the basis matrix has been successfully factorized */
      lpf->valid = 1;
#if _GLPLPF_DEBUG
      /* store the basis matrix for debugging */
      if (lpf->B != NULL) xfree(lpf->B);
      xassert(m <= 32767);
      lpf->B = B = xcalloc(1+m*m, sizeof(double));
      ind = xcalloc(1+m, sizeof(int));
      val = xcalloc(1+m, sizeof(double));
      for (k = 1; k <= m * m; k++)
         B[k] = 0.0;
      for (j = 1; j <= m; j++)
      {  len = col(info, j, ind, val);
         xassert(0 <= len && len <= m);
         for (k = 1; k <= len; k++)
         {  i = ind[k];
            xassert(1 <= i && i <= m);
            xassert(B[(i - 1) * m + j] == 0.0);
            xassert(val[k] != 0.0);
            B[(i - 1) * m + j] = val[k];
         }
      }
      xfree(ind);
      xfree(val);
#endif
      /* B = B0, so there are no additional rows/columns */
      lpf->n = 0;
      /* reset the Schur complement factorization */
      scf_reset_it(lpf->scf);
      /* P := Q := I */
      for (k = 1; k <= m; k++)
      {  lpf->P_row[k] = lpf->P_col[k] = k;
         lpf->Q_row[k] = lpf->Q_col[k] = k;
      }
      /* make all SVA locations free */
      lpf->v_ptr = 1;
      ret = 0;
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  The routine r_prod computes the product y := y + alpha * R * x,
*  where x is a n-vector, alpha is a scalar, y is a m0-vector.
*
*  Since matrix R is available by columns, the product is computed as
*  a linear combination:
*
*     y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]),
*
*  where R[j] is j-th column of R. */

static void r_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int j, beg, end, ptr;
      double t;
      for (j = 1; j <= n; j++)
      {  if (x[j] == 0.0) continue;
         /* y := y + alpha * R[j] * x[j] */
         t = a * x[j];
         beg = R_ptr[j];
         end = beg + R_len[j];
         for (ptr = beg; ptr < end; ptr++)
            y[v_ind[ptr]] += t * v_val[ptr];
      }
      return;
}

/***********************************************************************
*  The routine rt_prod computes the product y := y + alpha * R' * x,
*  where R' is a matrix transposed to R, x is a m0-vector, alpha is a
*  scalar, y is a n-vector.
*
*  Since matrix R is available by columns, the product components are
*  computed as inner products:
*
*     y[j] := y[j] + alpha * (j-th column of R) * x
*
*  for j = 1, 2, ..., n. */

static void rt_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int j, beg, end, ptr;
      double t;
      for (j = 1; j <= n; j++)
      {  /* t := (j-th column of R) * x */
         t = 0.0;
         beg = R_ptr[j];
         end = beg + R_len[j];
         for (ptr = beg; ptr < end; ptr++)
            t += v_val[ptr] * x[v_ind[ptr]];
         /* y[j] := y[j] + alpha * t */
         y[j] += a * t;
      }
      return;
}

/***********************************************************************
*  The routine s_prod computes the product y := y + alpha * S * x,
*  where x is a m0-vector, alpha is a scalar, y is a n-vector.
*
*  Since matrix S is available by rows, the product components are
*  computed as inner products:
*
*     y[i] = y[i] + alpha * (i-th row of S) * x
*
*  for i = 1, 2, ..., n. */

static void s_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *S_ptr = lpf->S_ptr;
      int *S_len = lpf->S_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int i, beg, end, ptr;
      double t;
      for (i = 1; i <= n; i++)
      {  /* t := (i-th row of S) * x */
         t = 0.0;
         beg = S_ptr[i];
         end = beg + S_len[i];
         for (ptr = beg; ptr < end; ptr++)
            t += v_val[ptr] * x[v_ind[ptr]];
         /* y[i] := y[i] + alpha * t */
         y[i] += a * t;
      }
      return;
}

/***********************************************************************
*  The routine st_prod computes the product y := y + alpha * S' * x,
*  where S' is a matrix transposed to S, x is a n-vector, alpha is a
*  scalar, y is m0-vector.
*
*  Since matrix R is available by rows, the product is computed as a
*  linear combination:
*
*     y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]),
*
*  where S'[i] is i-th row of S. */

static void st_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *S_ptr = lpf->S_ptr;
      int *S_len = lpf->S_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int i, beg, end, ptr;
      double t;
      for (i = 1; i <= n; i++)
      {  if (x[i] == 0.0) continue;
         /* y := y + alpha * S'[i] * x[i] */
         t = a * x[i];
         beg = S_ptr[i];
         end = beg + S_len[i];
         for (ptr = beg; ptr < end; ptr++)
            y[v_ind[ptr]] += t * v_val[ptr];
      }
      return;
}

#if _GLPLPF_DEBUG
/***********************************************************************
*  The routine check_error computes the maximal relative error between
*  left- and right-hand sides for the system B * x = b (if tr is zero)
*  or B' * x = b (if tr is non-zero), where B' is a matrix transposed
*  to B. (This routine is intended for debugging only.) */

static void check_error(LPF *lpf, int tr, const double x[],
      const double b[])
{     int m = lpf->m;
      double *B = lpf->B;
      int i, j;
      double  d, dmax = 0.0, s, t, tmax;
      for (i = 1; i <= m; i++)
      {  s = 0.0;
         tmax = 1.0;
         for (j = 1; j <= m; j++)
         {  if (!tr)
               t = B[m * (i - 1) + j] * x[j];
            else
               t = B[m * (j - 1) + i] * x[j];
            if (tmax < fabs(t)) tmax = fabs(t);
            s += t;
         }
         d = fabs(s - b[i]) / tmax;
         if (dmax < d) dmax = d;
      }
      if (dmax > 1e-8)
         xprintf("%s: dmax = %g; relative error too large\n",
            !tr ? "lpf_ftran" : "lpf_btran", dmax);
      return;
}
#endif

/***********************************************************************
*  NAME
*
*  lpf_ftran - perform forward transformation (solve system B*x = b)
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_ftran(LPF *lpf, double x[]);
*
*  DESCRIPTION
*
*  The routine lpf_ftran performs forward transformation, i.e. solves
*  the system B*x = b, where B is the basis matrix, x is the vector of
*  unknowns to be computed, b is the vector of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  BACKGROUND
*
*  Solution of the system B * x = b can be obtained by solving the
*  following augmented system:
*
*     ( B  F^) ( x )   ( b )
*     (      ) (   ) = (   )
*     ( G^ H^) ( y )   ( 0 )
*
*  which, using the main equality, can be written as follows:
*
*       ( L0 0 ) ( U0 R )   ( x )   ( b )
*     P (      ) (      ) Q (   ) = (   )
*       ( S  I ) ( 0  C )   ( y )   ( 0 )
*
*  therefore,
*
*     ( x )      ( U0 R )-1 ( L0 0 )-1    ( b )
*     (   ) = Q' (      )   (      )   P' (   )
*     ( y )      ( 0  C )   ( S  I )      ( 0 )
*
*  Thus, computing the solution includes the following steps:
*
*  1. Compute
*
*     ( f )      ( b )
*     (   ) = P' (   )
*     ( g )      ( 0 )
*
*  2. Solve the system
*
*     ( f1 )   ( L0 0 )-1 ( f )      ( L0 0 ) ( f1 )   ( f )
*     (    ) = (      )   (   )  =>  (      ) (    ) = (   )
*     ( g1 )   ( S  I )   ( g )      ( S  I ) ( g1 )   ( g )
*
*     from which it follows that:
*
*     { L0 * f1      = f      f1 = inv(L0) * f
*     {                   =>
*     {  S * f1 + g1 = g      g1 = g - S * f1
*
*  3. Solve the system
*
*     ( f2 )   ( U0 R )-1 ( f1 )      ( U0 R ) ( f2 )   ( f1 )
*     (    ) = (      )   (    )  =>  (      ) (    ) = (    )
*     ( g2 )   ( 0  C )   ( g1 )      ( 0  C ) ( g2 )   ( g1 )
*
*     from which it follows that:
*
*     { U0 * f2 + R * g2 = f1      f2 = inv(U0) * (f1 - R * g2)
*     {                        =>
*     {           C * g2 = g1      g2 = inv(C) * g1
*
*  4. Compute
*
*     ( x )      ( f2 )
*     (   ) = Q' (    )
*     ( y )      ( g2 )                                               */

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

/***********************************************************************
*  NAME
*
*  lpf_btran - perform backward transformation (solve system B'*x = b)
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_btran(LPF *lpf, double x[]);
*
*  DESCRIPTION
*
*  The routine lpf_btran performs backward transformation, i.e. solves
*  the system B'*x = b, where B' is a matrix transposed to the basis
*  matrix B, x is the vector of unknowns to be computed, b is the vector
*  of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  BACKGROUND
*
*  Solution of the system B' * x = b, where B' is a matrix transposed
*  to B, can be obtained by solving the following augmented system:
*
*     ( B  F^)T ( x )   ( b )
*     (      )  (   ) = (   )
*     ( G^ H^)  ( y )   ( 0 )
*
*  which, using the main equality, can be written as follows:
*
*      T ( U0 R )T ( L0 0 )T  T ( x )   ( b )
*     Q  (      )  (      )  P  (   ) = (   )
*        ( 0  C )  ( S  I )     ( y )   ( 0 )
*
*  or, equivalently, as follows:
*
*        ( U'0 0 ) ( L'0 S')    ( x )   ( b )
*     Q' (       ) (       ) P' (   ) = (   )
*        ( R'  C') ( 0   I )    ( y )   ( 0 )
*
*  therefore,
*
*     ( x )     ( L'0 S')-1 ( U'0 0 )-1   ( b )
*     (   ) = P (       )   (       )   Q (   )
*     ( y )     ( 0   I )   ( R'  C')     ( 0 )
*
*  Thus, computing the solution includes the following steps:
*
*  1. Compute
*
*     ( f )     ( b )
*     (   ) = Q (   )
*     ( g )     ( 0 )
*
*  2. Solve the system
*
*     ( f1 )   ( U'0 0 )-1 ( f )      ( U'0 0 ) ( f1 )   ( f )
*     (    ) = (       )   (   )  =>  (       ) (    ) = (   )
*     ( g1 )   ( R'  C')   ( g )      ( R'  C') ( g1 )   ( g )
*
*     from which it follows that:
*
*     { U'0 * f1           = f      f1 = inv(U'0) * f
*     {                         =>
*     { R'  * f1 + C' * g1 = g      g1 = inv(C') * (g - R' * f1)
*
*  3. Solve the system
*
*     ( f2 )   ( L'0 S')-1 ( f1 )      ( L'0 S') ( f2 )   ( f1 )
*     (    ) = (       )   (    )  =>  (       ) (    ) = (    )
*     ( g2 )   ( 0   I )   ( g1 )      ( 0   I ) ( g2 )   ( g1 )
*
*     from which it follows that:
*
*     { L'0 * f2 + S' * g2 = f1
*     {                          =>  f2 = inv(L'0) * ( f1 - S' * g2)
*     {                 g2 = g1
*
*  4. Compute
*
*     ( x )     ( f2 )
*     (   ) = P (    )
*     ( y )     ( g2 )                                                */

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 */
      luf_v_solve(lpf->luf, 1, f);
      /* 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);
      luf_f_solve(lpf->luf, 1, f);
      /* (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;
}

/***********************************************************************
*  The routine enlarge_sva enlarges the Sparse Vector Area to new_size
*  locations by reallocating the arrays v_ind and v_val. */

static void enlarge_sva(LPF *lpf, int new_size)
{     int v_size = lpf->v_size;
      int used = lpf->v_ptr - 1;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      xassert(v_size < new_size);
      while (v_size < new_size) v_size += v_size;
      lpf->v_size = v_size;
      lpf->v_ind = xcalloc(1+v_size, sizeof(int));
      lpf->v_val = xcalloc(1+v_size, sizeof(double));
      xassert(used >= 0);
      memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int));
      memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double));
      xfree(v_ind);
      xfree(v_val);
      return;
}

/***********************************************************************
*  NAME
*
*  lpf_update_it - update LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
*     const double val[]);
*
*  DESCRIPTION
*
*  The routine lpf_update_it updates the factorization of the basis
*  matrix B after replacing its j-th column by a new vector.
*
*  The parameter j specifies the number of column of B, which has been
*  replaced, 1 <= j <= m, where m is the order of B.
*
*  The parameter bh specifies the basis header entry for the new column
*  of B, which is the number of the new column in some original matrix.
*  This parameter is optional and can be specified as 0.
*
*  Row indices and numerical values of non-zero elements of the new
*  column of B should be placed in locations ind[1], ..., ind[len] and
*  val[1], ..., val[len], resp., where len is the number of non-zeros
*  in the column. Neither zero nor duplicate elements are allowed.
*
*  RETURNS
*
*  0  The factorization has been successfully updated.
*
*  LPF_ESING
*     New basis B is singular within the working precision.
*
*  LPF_ELIMIT
*     Maximal number of additional rows and columns has been reached.
*
*  BACKGROUND
*
*  Let j-th column of the current basis matrix B have to be replaced by
*  a new column a. This replacement is equivalent to removing the old
*  j-th column by fixing it at zero and introducing the new column as
*  follows:
*
*                   ( B   F^| a )
*     ( B  F^)      (       |   )
*     (      ) ---> ( G^  H^| 0 )
*     ( G^ H^)      (-------+---)
*                   ( e'j 0 | 0 )
*
*  where ej is a unit vector with 1 in j-th position which used to fix
*  the old j-th column of B (at zero). Then using the main equality we
*  have:
*
*     ( B   F^| a )            ( B0  F | f )
*     (       |   )   ( P  0 ) (       |   ) ( Q  0 )
*     ( G^  H^| 0 ) = (      ) ( G   H | g ) (      ) =
*     (-------+---)   ( 0  1 ) (-------+---) ( 0  1 )
*     ( e'j 0 | 0 )            ( v'  w'| 0 )
*
*       [   ( B0  F )|   ( f ) ]            [   ( B0 F )  |   ( f ) ]
*       [ P (       )| P (   ) ] ( Q  0 )   [ P (      ) Q| P (   ) ]
*     = [   ( G   H )|   ( g ) ] (      ) = [   ( G  H )  |   ( g ) ]
*       [------------+-------- ] ( 0  1 )   [-------------+---------]
*       [   ( v'  w')|     0   ]            [   ( v' w') Q|    0    ]
*
*  where:
*
*     ( a )     ( f )      ( f )        ( a )
*     (   ) = P (   )  =>  (   ) = P' * (   )
*     ( 0 )     ( g )      ( g )        ( 0 )
*
*                                 ( ej )      ( v )    ( v )     ( ej )
*     ( e'j  0 ) = ( v' w' ) Q => (    ) = Q' (   ) => (   ) = Q (    )
*                                 ( 0  )      ( w )    ( w )     ( 0  )
*
*  On the other hand:
*
*              ( B0| F  f )
*     ( P  0 ) (---+------) ( Q  0 )         ( B0    new F )
*     (      ) ( G | H  g ) (      ) = new P (             ) new Q
*     ( 0  1 ) (   |      ) ( 0  1 )         ( new G new H )
*              ( v'| w' 0 )
*
*  where:
*                               ( G )           ( H  g )
*     new F = ( F  f ), new G = (   ),  new H = (      ),
*                               ( v')           ( w' 0 )
*
*             ( P  0 )           ( Q  0 )
*     new P = (      ) , new Q = (      ) .
*             ( 0  1 )           ( 0  1 )
*
*  The factorization structure for the new augmented matrix remains the
*  same, therefore:
*
*           ( B0    new F )         ( L0     0 ) ( U0 new R )
*     new P (             ) new Q = (          ) (          )
*           ( new G new H )         ( new S  I ) ( 0  new C )
*
*  where:
*
*     new F = L0 * new R  =>
*
*     new R = inv(L0) * new F = inv(L0) * (F  f) = ( R  inv(L0)*f )
*
*     new G = new S * U0  =>
*
*                               ( G )             (     S      )
*     new S = new G * inv(U0) = (   ) * inv(U0) = (            )
*                               ( v')             ( v'*inv(U0) )
*
*     new H = new S * new R + new C  =>
*
*     new C = new H - new S * new R =
*
*             ( H  g )   (     S      )
*           = (      ) - (            ) * ( R  inv(L0)*f ) =
*             ( w' 0 )   ( v'*inv(U0) )
*
*             ( H - S*R           g - S*inv(L0)*f      )   ( C  x )
*           = (                                        ) = (      )
*             ( w'- v'*inv(U0)*R  -v'*inv(U0)*inv(L0)*f)   ( y' z )
*
*  Note that new C is resulted by expanding old C with new column x,
*  row y', and diagonal element z, where:
*
*     x = g - S * inv(L0) * f = g - S * (new column of R)
*
*     y = w - R'* inv(U'0)* v = w - R'* (new row of S)
*
*     z = - (new row of S) * (new column of R)
*
*  Finally, to replace old B by new B we have to permute j-th and last
*  (just added) columns of the matrix
*
*     ( B   F^| a )
*     (       |   )
*     ( G^  H^| 0 )
*     (-------+---)
*     ( e'j 0 | 0 )
*
*  and to keep the main equality do the same for matrix Q. */

int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
      const double val[])
{     int m0 = lpf->m0;
      int m = lpf->m;
#if _GLPLPF_DEBUG
      double *B = lpf->B;
#endif
      int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *S_ptr = lpf->S_ptr;
      int *S_len = lpf->S_len;
      int *P_row = lpf->P_row;
      int *P_col = lpf->P_col;
      int *Q_row = lpf->Q_row;
      int *Q_col = lpf->Q_col;
      int v_ptr = lpf->v_ptr;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      double *a = lpf->work2; /* new column */
      double *fg = lpf->work1, *f = fg, *g = fg + m0;
      double *vw = lpf->work2, *v = vw, *w = vw + m0;
      double *x = g, *y = w, z;
      int i, ii, k, ret;
      xassert(bh == bh);
      if (!lpf->valid)
         xfault("lpf_update_it: the factorization is not valid\n");
      if (!(1 <= j && j <= m))
         xfault("lpf_update_it: j = %d; column number out of range\n",
            j);
      xassert(0 <= m && m <= m0 + n);
      /* check if the basis factorization can be expanded */
      if (n == lpf->n_max)
      {  lpf->valid = 0;
         ret = LPF_ELIMIT;
         goto done;
      }
      /* convert new j-th column of B to dense format */
      for (i = 1; i <= m; i++)
         a[i] = 0.0;
      for (k = 1; k <= len; k++)
      {  i = ind[k];
         if (!(1 <= i && i <= m))
            xfault("lpf_update_it: ind[%d] = %d; row number out of rang"
               "e\n", k, i);
         if (a[i] != 0.0)
            xfault("lpf_update_it: ind[%d] = %d; duplicate row index no"
               "t allowed\n", k, i);
         if (val[k] == 0.0)
            xfault("lpf_update_it: val[%d] = %g; zero element not allow"
               "ed\n", k, val[k]);
         a[i] = val[k];
      }
#if _GLPLPF_DEBUG
      /* change column in the basis matrix for debugging */
      for (i = 1; i <= m; i++)
         B[(i - 1) * m + j] = a[i];
#endif
      /* (f g) := inv(P) * (a 0) */
      for (i = 1; i <= m0+n; i++)
         fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0);
      /* (v w) := Q * (ej 0) */
      for (i = 1; i <= m0+n; i++) vw[i] = 0.0;
      vw[Q_col[j]] = 1.0;
      /* f1 := inv(L0) * f (new column of R) */
      luf_f_solve(lpf->luf, 0, f);
      /* v1 := inv(U'0) * v (new row of S) */
      luf_v_solve(lpf->luf, 1, v);
      /* we need at most 2 * m0 available locations in the SVA to store
         new column of matrix R and new row of matrix S */
      if (lpf->v_size < v_ptr + m0 + m0)
      {  enlarge_sva(lpf, v_ptr + m0 + m0);
         v_ind = lpf->v_ind;
         v_val = lpf->v_val;
      }
      /* store new column of R */
      R_ptr[n+1] = v_ptr;
      for (i = 1; i <= m0; i++)
      {  if (f[i] != 0.0)
            v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++;
      }
      R_len[n+1] = v_ptr - lpf->v_ptr;
      lpf->v_ptr = v_ptr;
      /* store new row of S */
      S_ptr[n+1] = v_ptr;
      for (i = 1; i <= m0; i++)
      {  if (v[i] != 0.0)
            v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++;
      }
      S_len[n+1] = v_ptr - lpf->v_ptr;
      lpf->v_ptr = v_ptr;
      /* x := g - S * f1 (new column of C) */
      s_prod(lpf, x, -1.0, f);
      /* y := w - R' * v1 (new row of C) */
      rt_prod(lpf, y, -1.0, v);
      /* z := - v1 * f1 (new diagonal element of C) */
      z = 0.0;
      for (i = 1; i <= m0; i++) z -= v[i] * f[i];
      /* update factorization of new matrix C */
      switch (scf_update_exp(lpf->scf, x, y, z))
      {  case 0:
            break;
         case SCF_ESING:
            lpf->valid = 0;
            ret = LPF_ESING;
            goto done;
         case SCF_ELIMIT:
            xassert(lpf != lpf);
         default:
            xassert(lpf != lpf);
      }
      /* expand matrix P */
      P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1;
      /* expand matrix Q */
      Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1;
      /* permute j-th and last (just added) column of matrix Q */
      i = Q_col[j], ii = Q_col[m0+n+1];
      Q_row[i] = m0+n+1, Q_col[m0+n+1] = i;
      Q_row[ii] = j, Q_col[j] = ii;
      /* increase the number of additional rows and columns */
      lpf->n++;
      xassert(lpf->n <= lpf->n_max);
      /* the factorization has been successfully updated */
      ret = 0;
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  lpf_delete_it - delete LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_delete_it(LPF *lpf)
*
*  DESCRIPTION
*
*  The routine lpf_delete_it deletes LP basis factorization specified
*  by the parameter lpf and frees all memory allocated to this program
*  object. */

void lpf_delete_it(LPF *lpf)
{     luf_delete_it(lpf->luf);
#if _GLPLPF_DEBUG
      if (lpf->B != NULL) xfree(lpf->B);
#else
      xassert(lpf->B == NULL);
#endif
      if (lpf->R_ptr != NULL) xfree(lpf->R_ptr);
      if (lpf->R_len != NULL) xfree(lpf->R_len);
      if (lpf->S_ptr != NULL) xfree(lpf->S_ptr);
      if (lpf->S_len != NULL) xfree(lpf->S_len);
      if (lpf->scf != NULL) scf_delete_it(lpf->scf);
      if (lpf->P_row != NULL) xfree(lpf->P_row);
      if (lpf->P_col != NULL) xfree(lpf->P_col);
      if (lpf->Q_row != NULL) xfree(lpf->Q_row);
      if (lpf->Q_col != NULL) xfree(lpf->Q_col);
      if (lpf->v_ind != NULL) xfree(lpf->v_ind);
      if (lpf->v_val != NULL) xfree(lpf->v_val);
      if (lpf->work1 != NULL) xfree(lpf->work1);
      if (lpf->work2 != NULL) xfree(lpf->work2);
      xfree(lpf);
      return;
}

/* eof */