glpluf.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* glpluf.c (LU-factorization) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  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 "glplib.h"
#include "glpluf.h"
#define xfault xerror

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

#define N_MAX 100000000 /* = 100*10^6 */
/* maximal order of the original matrix */

/***********************************************************************
*  NAME
*
*  luf_create_it - create LU-factorization
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  LUF *luf_create_it(void);
*
*  DESCRIPTION
*
*  The routine luf_create_it creates a program object, which represents
*  LU-factorization of a square matrix.
*
*  RETURNS
*
*  The routine luf_create_it returns a pointer to the object created. */

LUF *luf_create_it(void)
{     LUF *luf;
      luf = xmalloc(sizeof(LUF));
      luf->n_max = luf->n = 0;
      luf->valid = 0;
      luf->fr_ptr = luf->fr_len = NULL;
      luf->fc_ptr = luf->fc_len = NULL;
      luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL;
      luf->vr_piv = NULL;
      luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL;
      luf->pp_row = luf->pp_col = NULL;
      luf->qq_row = luf->qq_col = NULL;
      luf->sv_size = 0;
      luf->sv_beg = luf->sv_end = 0;
      luf->sv_ind = NULL;
      luf->sv_val = NULL;
      luf->sv_head = luf->sv_tail = 0;
      luf->sv_prev = luf->sv_next = NULL;
      luf->vr_max = NULL;
      luf->rs_head = luf->rs_prev = luf->rs_next = NULL;
      luf->cs_head = luf->cs_prev = luf->cs_next = NULL;
      luf->flag = NULL;
      luf->work = NULL;
      luf->new_sva = 0;
      luf->piv_tol = 0.10;
      luf->piv_lim = 4;
      luf->suhl = 1;
      luf->eps_tol = 1e-15;
      luf->max_gro = 1e+10;
      luf->nnz_a = luf->nnz_f = luf->nnz_v = 0;
      luf->max_a = luf->big_v = 0.0;
      luf->rank = 0;
      return luf;
}

/***********************************************************************
*  NAME
*
*  luf_defrag_sva - defragment the sparse vector area
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_defrag_sva(LUF *luf);
*
*  DESCRIPTION
*
*  The routine luf_defrag_sva defragments the sparse vector area (SVA)
*  gathering all unused locations in one continuous extent. In order to
*  do that the routine moves all unused locations from the left part of
*  SVA (which contains rows and columns of the matrix V) to the middle
*  part (which contains free locations). This is attained by relocating
*  elements of rows and columns of the matrix V toward the beginning of
*  the left part.
*
*  NOTE that this "garbage collection" involves changing row and column
*  pointers of the matrix V. */

void luf_defrag_sva(LUF *luf)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_next = luf->sv_next;
      int sv_beg = 1;
      int i, j, k;
      /* skip rows and columns, which do not need to be relocated */
      for (k = luf->sv_head; k != 0; k = sv_next[k])
      {  if (k <= n)
         {  /* i-th row of the matrix V */
            i = k;
            if (vr_ptr[i] != sv_beg) break;
            vr_cap[i] = vr_len[i];
            sv_beg += vr_cap[i];
         }
         else
         {  /* j-th column of the matrix V */
            j = k - n;
            if (vc_ptr[j] != sv_beg) break;
            vc_cap[j] = vc_len[j];
            sv_beg += vc_cap[j];
         }
      }
      /* relocate other rows and columns in order to gather all unused
         locations in one continuous extent */
      for (k = k; k != 0; k = sv_next[k])
      {  if (k <= n)
         {  /* i-th row of the matrix V */
            i = k;
            memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]],
               vr_len[i] * sizeof(int));
            memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]],
               vr_len[i] * sizeof(double));
            vr_ptr[i] = sv_beg;
            vr_cap[i] = vr_len[i];
            sv_beg += vr_cap[i];
         }
         else
         {  /* j-th column of the matrix V */
            j = k - n;
            memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]],
               vc_len[j] * sizeof(int));
            memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]],
               vc_len[j] * sizeof(double));
            vc_ptr[j] = sv_beg;
            vc_cap[j] = vc_len[j];
            sv_beg += vc_cap[j];
         }
      }
      /* set new pointer to the beginning of the free part */
      luf->sv_beg = sv_beg;
      return;
}

/***********************************************************************
*  NAME
*
*  luf_enlarge_row - enlarge row capacity
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_enlarge_row(LUF *luf, int i, int cap);
*
*  DESCRIPTION
*
*  The routine luf_enlarge_row enlarges capacity of the i-th row of the
*  matrix V to cap locations (assuming that its current capacity is less
*  than cap). In order to do that the routine relocates elements of the
*  i-th row to the end of the left part of SVA (which contains rows and
*  columns of the matrix V) and then expands the left part by allocating
*  cap free locations from the free part. If there are less than cap
*  free locations, the routine defragments the sparse vector area.
*
*  Due to "garbage collection" this operation may change row and column
*  pointers of the matrix V.
*
*  RETURNS
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

int luf_enlarge_row(LUF *luf, int i, int cap)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      int ret = 0;
      int cur, k, kk;
      xassert(1 <= i && i <= n);
      xassert(vr_cap[i] < cap);
      /* if there are less than cap free locations, defragment SVA */
      if (luf->sv_end - luf->sv_beg < cap)
      {  luf_defrag_sva(luf);
         if (luf->sv_end - luf->sv_beg < cap)
         {  ret = 1;
            goto done;
         }
      }
      /* save current capacity of the i-th row */
      cur = vr_cap[i];
      /* copy existing elements to the beginning of the free part */
      memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]],
         vr_len[i] * sizeof(int));
      memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]],
         vr_len[i] * sizeof(double));
      /* set new pointer and new capacity of the i-th row */
      vr_ptr[i] = luf->sv_beg;
      vr_cap[i] = cap;
      /* set new pointer to the beginning of the free part */
      luf->sv_beg += cap;
      /* now the i-th row starts in the rightmost location among other
         rows and columns of the matrix V, so its node should be moved
         to the end of the row/column linked list */
      k = i;
      /* remove the i-th row node from the linked list */
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
      {  /* capacity of the previous row/column can be increased at the
            expense of old locations of the i-th row */
         kk = sv_prev[k];
         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
         sv_next[sv_prev[k]] = sv_next[k];
      }
      if (sv_next[k] == 0)
         luf->sv_tail = sv_prev[k];
      else
         sv_prev[sv_next[k]] = sv_prev[k];
      /* insert the i-th row node to the end of the linked list */
      sv_prev[k] = luf->sv_tail;
      sv_next[k] = 0;
      if (sv_prev[k] == 0)
         luf->sv_head = k;
      else
         sv_next[sv_prev[k]] = k;
      luf->sv_tail = k;
done: return ret;
}

/***********************************************************************
*  NAME
*
*  luf_enlarge_col - enlarge column capacity
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_enlarge_col(LUF *luf, int j, int cap);
*
*  DESCRIPTION
*
*  The routine luf_enlarge_col enlarges capacity of the j-th column of
*  the matrix V to cap locations (assuming that its current capacity is
*  less than cap). In order to do that the routine relocates elements
*  of the j-th column to the end of the left part of SVA (which contains
*  rows and columns of the matrix V) and then expands the left part by
*  allocating cap free locations from the free part. If there are less
*  than cap free locations, the routine defragments the sparse vector
*  area.
*
*  Due to "garbage collection" this operation may change row and column
*  pointers of the matrix V.
*
*  RETURNS
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

int luf_enlarge_col(LUF *luf, int j, int cap)
{     int n = luf->n;
      int *vr_cap = luf->vr_cap;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      int ret = 0;
      int cur, k, kk;
      xassert(1 <= j && j <= n);
      xassert(vc_cap[j] < cap);
      /* if there are less than cap free locations, defragment SVA */
      if (luf->sv_end - luf->sv_beg < cap)
      {  luf_defrag_sva(luf);
         if (luf->sv_end - luf->sv_beg < cap)
         {  ret = 1;
            goto done;
         }
      }
      /* save current capacity of the j-th column */
      cur = vc_cap[j];
      /* copy existing elements to the beginning of the free part */
      memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]],
         vc_len[j] * sizeof(int));
      memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]],
         vc_len[j] * sizeof(double));
      /* set new pointer and new capacity of the j-th column */
      vc_ptr[j] = luf->sv_beg;
      vc_cap[j] = cap;
      /* set new pointer to the beginning of the free part */
      luf->sv_beg += cap;
      /* now the j-th column starts in the rightmost location among
         other rows and columns of the matrix V, so its node should be
         moved to the end of the row/column linked list */
      k = n + j;
      /* remove the j-th column node from the linked list */
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
      {  /* capacity of the previous row/column can be increased at the