path: root/test/monniaux/glpk-4.65/src/bflib/fhv.c

/* fhv.c (sparse updatable FHV-factorization) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2012-2013 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 "env.h"
#include "fhv.h"

/***********************************************************************
*  fhv_ft_update - update FHV-factorization (Forrest-Tomlin)
*
*  This routine updates FHV-factorization of the original matrix A
*  after replacing its j-th column by a new one. The routine is based
*  on the method proposed by Forrest and Tomlin [1].
*
*  The parameter q specifies the number of column of A, which has been
*  replaced, 1 <= q <= n, where n is the order of A.
*
*  Row indices and numerical values of non-zero elements of the new
*  j-th column of A should be placed in locations aq_ind[1], ...,
*  aq_ind[aq_len] and aq_val[1], ..., aq_val[aq_len], respectively,
*  where aq_len is the number of non-zeros. Neither zero nor duplicate
*  elements are allowed.
*
*  The working arrays ind, val, and work should have at least 1+n
*  elements (0-th elements are not used).
*
*  RETURNS
*
*  0  The factorization has been successfully updated.
*
*  1  New matrix U = P'* V * Q' is upper triangular with zero diagonal
*     element u[s,s]. (Elimination was not performed.)
*
*  2  New matrix U = P'* V * Q' is upper triangular, and its diagonal
*     element u[s,s] or u[t,t] is too small in magnitude. (Elimination
*     was not performed.)
*
*  3  The same as 2, but after performing elimination.
*
*  4  The factorization has not been updated, because maximal number of
*     updates has been reached.
*
*  5  Accuracy test failed for the updated factorization.
*
*  BACKGROUND
*
*  The routine is based on the updating method proposed by Forrest and
*  Tomlin [1].
*
*  Let q-th column of the original matrix A have been replaced by new
*  column A[q]. Then, to keep the equality A = F * H * V, q-th column
*  of matrix V should be replaced by column V[q] = inv(F * H) * A[q].
*  From the standpoint of matrix U = P'* V * Q' such replacement is
*  equivalent to replacement of s-th column of matrix U, where s is
*  determined from q by permutation matrix Q. Thus, matrix U loses its
*  upper triangular form and becomes the following:
*
*        1   s       t   n
*     1  x x * x x x x x x
*        . x * x x x x x x
*     s  . . * x x x x x x
*        . . * x x x x x x
*        . . * . x x x x x
*        . . * . . x x x x
*     t  . . * . . . x x x
*        . . . . . . . x x
*     n  . . . . . . . . x
*
*  where t is largest row index of a non-zero element in s-th column.
*
*  The routine makes matrix U upper triangular as follows. First, it
*  moves rows and columns s+1, ..., t by one position to the left and
*  upwards, resp., and moves s-th row and s-th column to position t.
*  Due to such symmetric permutations matrix U becomes the following
*  (note that all diagonal elements remain on the diagonal, and element
*  u[s,s] becomes u[t,t]):
*
*        1   s       t   n
*     1  x x x x x x * x x
*        . x x x x x * x x
*     s  . . x x x x * x x
*        . . . x x x * x x
*        . . . . x x * x x
*        . . . . . x * x x
*     t  . . x x x x * x x
*        . . . . . . . x x
*     n  . . . . . . . . x
*
*  Then the routine performs gaussian elimination to eliminate
*  subdiagonal elements u[t,s], ..., u[t,t-1] using diagonal elements
*  u[s,s], ..., u[t-1,t-1] as pivots. During the elimination process
*  the routine permutes neither rows nor columns, so only t-th row is
*  changed. Should note that actually all operations are performed on
*  matrix V = P * U * Q, since matrix U is not stored.
*
*  To keep the equality A = F * H * V, the routine appends new row-like
*  factor H[k] to matrix H, and every time it applies elementary
*  gaussian transformation to eliminate u[t,j'] = v[p,j] using pivot
*  u[j',j'] = v[i,j], it also adds new element f[p,j] = v[p,j] / v[i,j]
*  (gaussian multiplier) to factor H[k], which initially is a unity
*  matrix. At the end of elimination process the row-like factor H[k]
*  may look as follows:
*
*        1               n          1   s       t   n
*     1  1 . . . . . . . .       1  1 . . . . . . . .
*        . 1 . . . . . . .          . 1 . . . . . . .
*        . . 1 . . . . . .       s  . . 1 . . . . . .
*     p  . x x 1 . x . x .          . . . 1 . . . . .
*        . . . . 1 . . . .          . . . . 1 . . . .
*        . . . . . 1 . . .          . . . . . 1 . . .
*        . . . . . . 1 . .       t  . . x x x x 1 . .
*        . . . . . . . 1 .          . . . . . . . 1 .
*     n  . . . . . . . . 1       n  . . . . . . . . 1
*
*              H[k]                 inv(P) * H[k] * P
*
*  If, however, s = t, no elimination is needed, in which case no new
*  row-like factor is created.
*
*  REFERENCES
*
*  1. J.J.H.Forrest and J.A.Tomlin, "Updated triangular factors of the
*     basis to maintain sparsity in the product form simplex method,"
*     Math. Prog. 2 (1972), pp. 263-78. */

int fhv_ft_update(FHV *fhv, int q, int aq_len, const int aq_ind[],
      const double aq_val[], int ind[/*1+n*/], double val[/*1+n*/],
      double work[/*1+n*/])
{     LUF *luf = fhv->luf;
      int n = luf->n;
      SVA *sva = luf->sva;
      int *sv_ind = sva->ind;
      double *sv_val = sva->val;
      int vr_ref = luf->vr_ref;
      int *vr_ptr = &sva->ptr[vr_ref-1];
      int *vr_len = &sva->len[vr_ref-1];
      int *vr_cap = &sva->cap[vr_ref-1];
      double *vr_piv = luf->vr_piv;
      int vc_ref = luf->vc_ref;
      int *vc_ptr = &sva->ptr[vc_ref-1];
      int *vc_len = &sva->len[vc_ref-1];
      int *vc_cap = &sva->cap[vc_ref-1];
      int *pp_ind = luf->pp_ind;
      int *pp_inv = luf->pp_inv;
      int *qq_ind = luf->qq_ind;
      int *qq_inv = luf->qq_inv;
      int *hh_ind = fhv->hh_ind;
      int hh_ref = fhv->hh_ref;
      int *hh_ptr = &sva->ptr[hh_ref-1];
      int *hh_len = &sva->len[hh_ref-1];
#if 1 /* FIXME */
      const double eps_tol = DBL_EPSILON;
      const double vpq_tol = 1e-5;
      const double err_tol = 1e-10;
#endif
      int end, i, i_end, i_ptr, j, j_end, j_ptr, k, len, nnz, p, p_end,
         p_ptr, ptr, q_end, q_ptr, s, t;
      double f, vpq, temp;
      /*--------------------------------------------------------------*/
      /* replace current q-th column of matrix V by new one           */
      /*--------------------------------------------------------------*/
      xassert(1 <= q && q <= n);
      /* convert new q-th column of matrix A to dense format */
      for (i = 1; i <= n; i++)
         val[i] = 0.0;
      xassert(0 <= aq_len && aq_len <= n);
      for (k = 1; k <= aq_len; k++)
      {  i = aq_ind[k];
         xassert(1 <= i && i <= n);
         xassert(val[i] == 0.0);
         xassert(aq_val[k] != 0.0);
         val[i] = aq_val[k];
      }
      /* compute new q-th column of matrix V:
       * new V[q] = inv(F * H) * (new A[q]) */
      luf->pp_ind = fhv->p0_ind;
      luf->pp_inv = fhv->p0_inv;
      luf_f_solve(luf, val);
      luf->pp_ind = pp_ind;
      luf->pp_inv = pp_inv;
      fhv_h_solve(fhv, val);
      /* q-th column of V = s-th column of U */
      s = qq_inv[q];
      /* determine row number of element v[p,q] that corresponds to
       * diagonal element u[s,s] */
      p = pp_inv[s];
      /* convert new q-th column of V to sparse format;
       * element v[p,q] = u[s,s] is not included in the element list
       * and stored separately */
      vpq = 0.0;
      len = 0;
      for (i = 1; i <= n; i++)
      {  temp = val[i];
#if 1 /* FIXME */
         if (-eps_tol < temp && temp < +eps_tol)
#endif
            /* nop */;
         else if (i == p)
            vpq = temp;
         else
         {  ind[++len] = i;
            val[len] = temp;
         }
      }
      /* clear q-th column of matrix V */
      for (q_end = (q_ptr = vc_ptr[q]) + vc_len[q];
         q_ptr < q_end; q_ptr++)
      {  /* get row index of v[i,q] */
         i = sv_ind[q_ptr];
         /* find and remove v[i,q] from i-th row */
         for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i];
            sv_ind[i_ptr] != q; i_ptr++)
            /* nop */;
         xassert(i_ptr < i_end);
         sv_ind[i_ptr] = sv_ind[i_end-1];
         sv_val[i_ptr] = sv_val[i_end-1];
         vr_len[i]--;
      }
      /* now q-th column of matrix V is empty */
      vc_len[q] = 0;
      /* put new q-th column of V (except element v[p,q] = u[s,s]) in
       * column-wise format */
      if (len > 0)
      {  if (vc_cap[q] < len)
         {  if (sva->r_ptr - sva->m_ptr < len)
            {  sva_more_space(sva, len);
               sv_ind = sva->ind;
               sv_val = sva->val;
            }
            sva_enlarge_cap(sva, vc_ref-1+q, len, 0);
         }
         ptr = vc_ptr[q];
         memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int));
         memcpy(&sv_val[ptr], &val[1], len * sizeof(double));
         vc_len[q] = len;
      }
      /* put new q-th column of V (except element v[p,q] = u[s,s]) in
       * row-wise format, and determine largest row number t such that
       * u[s,t] != 0 */
      t = (vpq == 0.0 ? 0 : s);
      for (k = 1; k <= len; k++)
      {  /* get row index of v[i,q] */
         i = ind[k];
         /* put v[i,q] to i-th row */
         if (vr_cap[i] == vr_len[i])
         {  /* reserve extra locations in i-th row to reduce further
             * relocations of that row */
#if 1 /* FIXME */
            int need = vr_len[i] + 5;
#endif
            if (sva->r_ptr - sva->m_ptr < need)
            {  sva_more_space(sva, need);
               sv_ind = sva->ind;
               sv_val = sva->val;
            }
            sva_enlarge_cap(sva, vr_ref-1+i, need, 0);
         }
         sv_ind[ptr = vr_ptr[i] + (vr_len[i]++)] = q;
         sv_val[ptr] = val[k];
         /* v[i,q] is non-zero; increase t */
         if (t < pp_ind[i])
            t = pp_ind[i];
      }
      /*--------------------------------------------------------------*/
      /* check if matrix U is already upper triangular                */
      /*--------------------------------------------------------------*/
      /* check if there is a spike in s-th column of matrix U, which
       * is q-th column of matrix V */
      if (s >= t)
      {  /* no spike; matrix U is already upper triangular */
         /* store its diagonal element u[s,s] = v[p,q] */
         vr_piv[p] = vpq;
         if (s > t)
         {  /* matrix U is structurally singular, because its diagonal
             * element u[s,s] = v[p,q] is exact zero */
            xassert(vpq == 0.0);
            return 1;
         }
#if 1 /* FIXME */
         else if (-vpq_tol < vpq && vpq < +vpq_tol)
#endif
         {  /* matrix U is not well conditioned, because its diagonal
             * element u[s,s] = v[p,q] is too small in magnitude */
            return 2;
         }
         else
         {  /* normal case */
            return 0;
         }
      }
      /*--------------------------------------------------------------*/
      /* perform implicit symmetric permutations of rows and columns  */
      /* of matrix U                                                  */
      /*--------------------------------------------------------------*/
      /* currently v[p,q] = u[s,s] */
      xassert(p == pp_inv[s] && q == qq_ind[s]);
      for (k = s; k < t; k++)
      {  pp_ind[pp_inv[k] = pp_inv[k+1]] = k;
         qq_inv[qq_ind[k] = qq_ind[k+1]] = k;
      }
      /* now v[p,q] = u[t,t] */
      pp_ind[pp_inv[t] = p] = qq_inv[qq_ind[t] = q] = t;
      /*--------------------------------------------------------------*/
      /* check if matrix U is already upper triangular                */
      /*--------------------------------------------------------------*/
      /* check if there is a spike in t-th row of matrix U, which is
       * p-th row of matrix V */
      for (p_end = (p_ptr = vr_ptr[p]) + vr_len[p];
         p_ptr < p_end; p_ptr++)
      {  if (qq_inv[sv_ind[p_ptr]] < t)
            break; /* spike detected */
      }
      if (p_ptr == p_end)
      {  /* no spike; matrix U is already upper triangular */
         /* store its diagonal element u[t,t] = v[p,q] */
         vr_piv[p] = vpq;
#if 1 /* FIXME */
         if (-vpq_tol < vpq && vpq < +vpq_tol)
#endif
         {  /* matrix U is not well conditioned, because its diagonal
             * element u[t,t] = v[p,q] is too small in magnitude */
            return 2;
         }
         else
         {  /* normal case */
            return 0;
         }
      }
      /*--------------------------------------------------------------*/
      /* copy p-th row of matrix V, which is t-th row of matrix U, to */
      /* working array                                                */
      /*--------------------------------------------------------------*/
      /* copy p-th row of matrix V, including element v[p,q] = u[t,t],
       * to the working array in dense format and remove these elements
       * from matrix V; since no pivoting is used, only this row will
       * change during elimination */
      for (j = 1; j <= n; j++)
         work[j] = 0.0;
      work[q] = vpq;
      for (p_end = (p_ptr = vr_ptr[p]) + vr_len[p];
         p_ptr < p_end; p_ptr++)
      {  /* get column index of v[p,j] and store this element to the
          * working array */
         work[j = sv_ind[p_ptr]] = sv_val[p_ptr];
         /* find and remove v[p,j] from j-th column */
         for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j];
            sv_ind[j_ptr] != p; j_ptr++)
            /* nop */;
         xassert(j_ptr < j_end);
         sv_ind[j_ptr] = sv_ind[j_end-1];
         sv_val[j_ptr] = sv_val[j_end-1];
         vc_len[j]--;
      }
      /* now p-th row of matrix V is temporarily empty */
      vr_len[p] = 0;
      /*--------------------------------------------------------------*/
      /* perform gaussian elimination                                 */
      /*--------------------------------------------------------------*/
      /* transform p-th row of matrix V stored in working array, which
       * is t-th row of matrix U, to eliminate subdiagonal elements
       * u[t,s], ..., u[t,t-1]; corresponding gaussian multipliers will
       * form non-trivial row of new row-like factor */
      nnz = 0; /* number of non-zero gaussian multipliers */
      for (k = s; k < t; k++)
      {  /* diagonal element u[k,k] = v[i,j] is used as pivot */
         i = pp_inv[k], j = qq_ind[k];
         /* take subdiagonal element u[t,k] = v[p,j] */
         temp = work[j];
#if 1 /* FIXME */
         if (-eps_tol < temp && temp < +eps_tol)
            continue;
#endif
         /* compute and save gaussian multiplier:
          * f := u[t,k] / u[k,k] = v[p,j] / v[i,j] */
         ind[++nnz] = i;
         val[nnz] = f = work[j] / vr_piv[i];
         /* gaussian transformation to eliminate u[t,k] = v[p,j]:
          * (p-th row of V) := (p-th row of V) - f * (i-th row of V) */
         for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i];
            i_ptr < i_end; i_ptr++)
            work[sv_ind[i_ptr]] -= f * sv_val[i_ptr];
      }
      /* now matrix U is again upper triangular */
#if 1 /* FIXME */
      if (-vpq_tol < work[q] && work[q] < +vpq_tol)
#endif
      {  /* however, its new diagonal element u[t,t] = v[p,q] is too
          * small in magnitude */
         return 3;
      }
      /*--------------------------------------------------------------*/
      /* create new row-like factor H[k] and add to eta file H        */
      /*--------------------------------------------------------------*/
      /* (nnz = 0 means that all subdiagonal elements were too small
       * in magnitude) */
      if (nnz > 0)
      {  if (fhv->nfs == fhv->nfs_max)
         {  /* maximal number of row-like factors has been reached */
            return 4;
         }
         k = ++(fhv->nfs);
         hh_ind[k] = p;
         /* store non-trivial row of H[k] in right (dynamic) part of
          * SVA (diagonal unity element is not stored) */
         if (sva->r_ptr - sva->m_ptr < nnz)
         {  sva_more_space(sva, nnz);
            sv_ind = sva->ind;
            sv_val = sva->val;
         }
         sva_reserve_cap(sva, fhv->hh_ref-1+k, nnz);
         ptr = hh_ptr[k];
         memcpy(&sv_ind[ptr], &ind[1], nnz * sizeof(int));
         memcpy(&sv_val[ptr], &val[1], nnz * sizeof(double));
         hh_len[k] = nnz;
      }
      /*--------------------------------------------------------------*/
      /* copy transformed p-th row of matrix V, which is t-th row of  */
      /* matrix U, from working array back to matrix V                */
      /*--------------------------------------------------------------*/
      /* copy elements of transformed p-th row of matrix V, which are
       * non-diagonal elements u[t,t+1], ..., u[t,n] of matrix U, from
       * working array to corresponding columns of matrix V (note that
       * diagonal element u[t,t] = v[p,q] not copied); also transform
       * p-th row of matrix V to sparse format */
      len = 0;
      for (k = t+1; k <= n; k++)
      {  /* j-th column of V = k-th column of U */
         j = qq_ind[k];
         /* take non-diagonal element v[p,j] = u[t,k] */
         temp = work[j];
#if 1 /* FIXME */
         if (-eps_tol < temp && temp < +eps_tol)
            continue;
#endif
         /* add v[p,j] to j-th column of matrix V */
         if (vc_cap[j] == vc_len[j])
         {  /* reserve extra locations in j-th column to reduce further
             * relocations of that column */
#if 1 /* FIXME */
            int need = vc_len[j] + 5;
#endif
            if (sva->r_ptr - sva->m_ptr < need)
            {  sva_more_space(sva, need);
               sv_ind = sva->ind;
               sv_val = sva->val;
            }
            sva_enlarge_cap(sva, vc_ref-1+j, need, 0);
         }
         sv_ind[ptr = vc_ptr[j] + (vc_len[j]++)] = p;
         sv_val[ptr] = temp;
         /* store element v[p,j] = u[t,k] to working sparse vector */
         ind[++len] = j;
         val[len] = temp;
      }
      /* copy elements from working sparse vector to p-th row of matrix
       * V (this row is currently empty) */
      if (vr_cap[p] < len)
      {  if (sva->r_ptr - sva->m_ptr < len)
         {  sva_more_space(sva, len);
            sv_ind = sva->ind;
            sv_val = sva->val;
         }
         sva_enlarge_cap(sva, vr_ref-1+p, len, 0);
      }
      ptr = vr_ptr[p];
      memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int));
      memcpy(&sv_val[ptr], &val[1], len * sizeof(double));
      vr_len[p] = len;
      /* store new diagonal element u[t,t] = v[p,q] */
      vr_piv[p] = work[q];
      /*--------------------------------------------------------------*/
      /* perform accuracy test (only if new H[k] was added)           */
      /*--------------------------------------------------------------*/
      if (nnz > 0)
      {  /* copy p-th (non-trivial) row of row-like factor H[k] (except
          * unity diagonal element) to working array in dense format */
         for (j = 1; j <= n; j++)
            work[j] = 0.0;
         k = fhv->nfs;
         for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++)
            work[sv_ind[ptr]] = sv_val[ptr];
         /* compute inner product of p-th (non-trivial) row of matrix
          * H[k] and q-th column of matrix V */
         temp = vr_piv[p]; /* 1 * v[p,q] */
         ptr = vc_ptr[q];
         end = ptr + vc_len[q];
         for (; ptr < end; ptr++)
            temp += work[sv_ind[ptr]] * sv_val[ptr];
         /* inner product should be equal to element v[p,q] *before*
          * matrix V was transformed */
         /* compute relative error */
         temp = fabs(vpq - temp) / (1.0 + fabs(vpq));
#if 1 /* FIXME */
         if (temp > err_tol)
#endif
         {  /* relative error is too large */
            return 5;
         }
      }
      /* factorization has been successfully updated */
      return 0;
}

/***********************************************************************
*  fhv_h_solve - solve system H * x = b
*
*  This routine solves the system H * x = b, where the matrix H is the
*  middle factor of the sparse updatable FHV-factorization.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[n], where n is the order of the
*  matrix H. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void fhv_h_solve(FHV *fhv, double x[/*1+n*/])
{     SVA *sva = fhv->luf->sva;
      int *sv_ind = sva->ind;
      double *sv_val = sva->val;
      int nfs = fhv->nfs;
      int *hh_ind = fhv->hh_ind;
      int hh_ref = fhv->hh_ref;
      int *hh_ptr = &sva->ptr[hh_ref-1];
      int *hh_len = &sva->len[hh_ref-1];
      int i, k, end, ptr;
      double x_i;
      for (k = 1; k <= nfs; k++)
      {  x_i = x[i = hh_ind[k]];
         for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++)
            x_i -= sv_val[ptr] * x[sv_ind[ptr]];
         x[i] = x_i;
      }
      return;
}

/***********************************************************************
*  fhv_ht_solve - solve system H' * x = b
*
*  This routine solves the system H' * x = b, where H' is a matrix
*  transposed to the matrix H, which is the middle factor of the sparse
*  updatable FHV-factorization.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[n], where n is the order of the
*  matrix H. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void fhv_ht_solve(FHV *fhv, double x[/*1+n*/])
{     SVA *sva = fhv->luf->sva;
      int *sv_ind = sva->ind;
      double *sv_val = sva->val;
      int nfs = fhv->nfs;
      int *hh_ind = fhv->hh_ind;
      int hh_ref = fhv->hh_ref;
      int *hh_ptr = &sva->ptr[hh_ref-1];
      int *hh_len = &sva->len[hh_ref-1];
      int k, end, ptr;
      double x_j;
      for (k = nfs; k >= 1; k--)
      {  if ((x_j = x[hh_ind[k]]) == 0.0)
            continue;
         for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++)
            x[sv_ind[ptr]] -= sv_val[ptr] * x_j;
      }
      return;
}

/* eof */