/* 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 */