glpspm.c

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

*  RETURNS
*
*  The routine returns the number of elements removed. */

int spm_drop_zeros(SPM *A, double eps)
{     SPME *e, *next;
      int i, count = 0;
      for (i = 1; i <= A->m; i++)
      {  for (e = A->row[i]; e != NULL; e = next)
         {  next = e->r_next;
            if (e->val == 0.0 || fabs(e->val) < eps)
            {  /* remove element from the row list */
               if (e->r_prev == NULL)
                  A->row[e->i] = e->r_next;
               else
                  e->r_prev->r_next = e->r_next;
               if (e->r_next == NULL)
                  ;
               else
                  e->r_next->r_prev = e->r_prev;
               /* remove element from the column list */
               if (e->c_prev == NULL)
                  A->col[e->j] = e->c_next;
               else
                  e->c_prev->c_next = e->c_next;
               if (e->c_next == NULL)
                  ;
               else
                  e->c_next->c_prev = e->c_prev;
               /* return element to the memory pool */
               dmp_free_atom(A->pool, e, sizeof(SPME));
               count++;
            }
         }
      }
      return count;
}

/***********************************************************************
*  NAME
*
*  spm_read_mat - read sparse matrix from text file
*
*  SYNOPSIS
*
*  #include "glpspm.h"
*  SPM *spm_read_mat(const char *fname);
*
*  DESCRIPTION
*
*  The routine reads a sparse matrix from a text file whose name is
*  specified by the parameter fname.
*
*  For the file format see description of the routine spm_write_mat.
*
*  RETURNS
*
*  On success the routine returns a pointer to the matrix created,
*  otherwise NULL. */

SPM *spm_read_mat(const char *fname)
{     SPM *A = NULL;
      PDS *pds;
      jmp_buf jump;
      int i, j, k, m, n, nnz, fail = 0;
      double val;
      xprintf("spm_read_mat: reading matrix from `%s'...\n", fname);
      pds = pds_open_file(fname);
      if (pds == NULL)
      {  xprintf("spm_read_mat: unable to open `%s' - %s\n", fname,
            strerror(errno));
         fail = 1;
         goto done;
      }
      if (setjmp(jump))
      {  fail = 1;
         goto done;
      }
      pds_set_jump(pds, jump);
      /* number of rows, number of columns, number of non-zeros */
      m = pds_scan_int(pds);
      if (m < 0)
         pds_error(pds, "invalid number of rows\n");
      n = pds_scan_int(pds);
      if (n < 0)
         pds_error(pds, "invalid number of columns\n");
      nnz = pds_scan_int(pds);
      if (nnz < 0)
         pds_error(pds, "invalid number of non-zeros\n");
      /* create matrix */
      xprintf("spm_read_mat: %d rows, %d columns, %d non-zeros\n",
         m, n, nnz);
      A = spm_create_mat(m, n);
      /* read matrix elements */
      for (k = 1; k <= nnz; k++)
      {  /* row index, column index, element value */
         i = pds_scan_int(pds);
         if (!(1 <= i && i <= m))
            pds_error(pds, "row index out of range\n");
         j = pds_scan_int(pds);
         if (!(1 <= j && j <= n))
            pds_error(pds, "column index out of range\n");
         val = pds_scan_num(pds);
         /* add new element to the matrix */
         spm_new_elem(A, i, j, val);
      }
      xprintf("spm_read_mat: %d lines were read\n", pds->count);
done: if (pds != NULL) pds_close_file(pds);
      if (fail && A != NULL) spm_delete_mat(A), A = NULL;
      return A;
}

/***********************************************************************
*  NAME
*
*  spm_write_mat - write sparse matrix to text file
*
*  SYNOPSIS
*
*  #include "glpspm.h"
*  int spm_write_mat(const SPM *A, const char *fname);
*
*  DESCRIPTION
*
*  The routine spm_write_mat writes the specified sparse matrix to a
*  text file whose name is specified by the parameter fname. This file
*  can be read back with the routine spm_read_mat.
*
*  RETURNS
*
*  On success the routine returns zero, otherwise non-zero.
*
*  FILE FORMAT
*
*  The file created by the routine spm_write_mat is a plain text file,
*  which contains the following information:
*
*     m n nnz
*     row[1] col[1] val[1]
*     row[2] col[2] val[2]
*     . . .
*     row[nnz] col[nnz] val[nnz]
*
*  where:
*  m is the number of rows;
*  n is the number of columns;
*  nnz is the number of non-zeros;
*  row[k], k = 1,...,nnz, are row indices;
*  col[k], k = 1,...,nnz, are column indices;
*  val[k], k = 1,...,nnz, are element values. */

int spm_write_mat(const SPM *A, const char *fname)
{     FILE *fp;
      int i, nnz, ret = 0;
      xprintf("spm_write_mat: writing matrix to `%s'...\n", fname);
      fp = fopen(fname, "w");
      if (fp == NULL)
      {  xprintf("spm_write_mat: unable to create `%s' - %s\n", fname,
            strerror(errno));
         ret = 1;
         goto done;
      }
      /* number of rows, number of columns, number of non-zeros */
      nnz = spm_count_nnz(A);
      fprintf(fp, "%d %d %d\n", A->m, A->n, nnz);
      /* walk through rows of the matrix */
      for (i = 1; i <= A->m; i++)
      {  SPME *e;
         /* walk through elements of i-th row */
         for (e = A->row[i]; e != NULL; e = e->r_next)
         {  /* row index, column index, element value */
            fprintf(fp, "%d %d %.*g\n", e->i, e->j, DBL_DIG, e->val);
         }
      }
      fflush(fp);
      if (ferror(fp))
      {  xprintf("spm_write_mat: writing error on `%s' - %s\n", fname,
            strerror(errno));
         ret = 1;
         goto done;
      }
      xprintf("spm_write_mat: %d lines were written\n", 1 + nnz);
done: if (fp != NULL) fclose(fp);
      return ret;
}

/***********************************************************************
*  NAME
*
*  spm_transpose - transpose sparse matrix
*
*  SYNOPSIS
*
*  #include "glpspm.h"
*  SPM *spm_transpose(const SPM *A);
*
*  RETURNS
*
*  The routine computes and returns sparse matrix B, which is a matrix
*  transposed to sparse matrix A. */

SPM *spm_transpose(const SPM *A)
{     SPM *B;
      int i;
      B = spm_create_mat(A->n, A->m);
      for (i = 1; i <= A->m; i++)
      {  SPME *e;
         for (e = A->row[i]; e != NULL; e = e->r_next)
            spm_new_elem(B, e->j, i, e->val);
      }
      return B;
}

SPM *spm_add_sym(const SPM *A, const SPM *B)
{     /* add two sparse matrices (symbolic phase) */
      SPM *C;
      int i, j, *flag;
      xassert(A->m == B->m);
      xassert(A->n == B->n);
      /* create resultant matrix */
      C = spm_create_mat(A->m, A->n);
      /* allocate and clear the flag array */
      flag = xcalloc(1+C->n, sizeof(int));
      for (j = 1; j <= C->n; j++)
         flag[j] = 0;
      /* compute pattern of C = A + B */
      for (i = 1; i <= C->m; i++)
      {  SPME *e;
         /* at the beginning i-th row of C is empty */
         /* (i-th row of C) := (i-th row of C) union (i-th row of A) */
         for (e = A->row[i]; e != NULL; e = e->r_next)
         {  /* (note that i-th row of A may have duplicate elements) */
            j = e->j;
            if (!flag[j])
            {  spm_new_elem(C, i, j, 0.0);
               flag[j] = 1;
            }
         }
         /* (i-th row of C) := (i-th row of C) union (i-th row of B) */
         for (e = B->row[i]; e != NULL; e = e->r_next)
         {  /* (note that i-th row of B may have duplicate elements) */
            j = e->j;
            if (!flag[j])
            {  spm_new_elem(C, i, j, 0.0);
               flag[j] = 1;
            }
         }
         /* reset the flag array */
         for (e = C->row[i]; e != NULL; e = e->r_next)
            flag[e->j] = 0;
      }
      /* check and deallocate the flag array */
      for (j = 1; j <= C->n; j++)
         xassert(!flag[j]);
      xfree(flag);
      return C;
}

void spm_add_num(SPM *C, double alfa, const SPM *A, double beta,
      const SPM *B)
{     /* add two sparse matrices (numeric phase) */
      int i, j;
      double *work;
      /* allocate and clear the working array */
      work = xcalloc(1+C->n, sizeof(double));
      for (j = 1; j <= C->n; j++)
         work[j] = 0.0;
      /* compute matrix C = alfa * A + beta * B */
      for (i = 1; i <= C->n; i++)
      {  SPME *e;
         /* work := alfa * (i-th row of A) + beta * (i-th row of B) */
         /* (note that A and/or B may have duplicate elements) */
         for (e = A->row[i]; e != NULL; e = e->r_next)
            work[e->j] += alfa * e->val;
         for (e = B->row[i]; e != NULL; e = e->r_next)
            work[e->j] += beta * e->val;
         /* (i-th row of C) := work, work := 0 */
         for (e = C->row[i]; e != NULL; e = e->r_next)
         {  j = e->j;
            e->val = work[j];
            work[j] = 0.0;
         }
      }
      /* check and deallocate the working array */
      for (j = 1; j <= C->n; j++)
         xassert(work[j] == 0.0);
      xfree(work);
      return;
}

SPM *spm_add_mat(double alfa, const SPM *A, double beta, const SPM *B)
{     /* add two sparse matrices (driver routine) */
      SPM *C;
      C = spm_add_sym(A, B);
      spm_add_num(C, alfa, A, beta, B);
      return C;
}

SPM *spm_mul_sym(const SPM *A, const SPM *B)
{     /* multiply two sparse matrices (symbolic phase) */
      int i, j, k, *flag;
      SPM *C;
      xassert(A->n == B->m);
      /* create resultant matrix */
      C = spm_create_mat(A->m, B->n);
      /* allocate and clear the flag array */
      flag = xcalloc(1+C->n, sizeof(int));
      for (j = 1; j <= C->n; j++)
         flag[j] = 0;
      /* compute pattern of C = A * B */
      for (i = 1; i <= C->m; i++)
      {  SPME *e, *ee;
         /* compute pattern of i-th row of C */
         for (e = A->row[i]; e != NULL; e = e->r_next)
         {  k = e->j;
            for (ee = B->row[k]; ee != NULL; ee = ee->r_next)
            {  j = ee->j;
               /* if a[i,k] != 0 and b[k,j] != 0 then c[i,j] != 0 */
               if (!flag[j])
               {  /* c[i,j] does not exist, so create it */
                  spm_new_elem(C, i, j, 0.0);
                  flag[j] = 1;
               }
            }
         }
         /* reset the flag array */
         for (e = C->row[i]; e != NULL; e = e->r_next)
            flag[e->j] = 0;
      }
      /* check and deallocate the flag array */
      for (j = 1; j <= C->n; j++)
         xassert(!flag[j]);
      xfree(flag);
      return C;
}

void spm_mul_num(SPM *C, const SPM *A, const SPM *B)
{     /* multiply two sparse matrices (numeric phase) */
      int i, j;
      double *work;
      /* allocate and clear the working array */
      work = xcalloc(1+A->n, sizeof(double));
      for (j = 1; j <= A->n; j++)
         work[j] = 0.0;
      /* compute matrix C = A * B */
      for (i = 1; i <= C->m; i++)
      {  SPME *e, *ee;
         double temp;
         /* work := (i-th row of A) */
         /* (note that A may have duplicate elements) */
         for (e = A->row[i]; e != NULL; e = e->r_next)
            work[e->j] += e->val;
         /* compute i-th row of C */
         for (e = C->row[i]; e != NULL; e = e->r_next)
         {  j = e->j;
            /* c[i,j] := work * (j-th column of B) */
            temp = 0.0;
            for (ee = B->col[j]; ee != NULL; ee = ee->c_next)
               temp += work[ee->i] * ee->val;
            e->val = temp;
         }
         /* reset the working array */
         for (e = A->row[i]; e != NULL; e = e->r_next)
            work[e->j] = 0.0;
      }
      /* check and deallocate the working array */
      for (j = 1; j <= A->n; j++)
         xassert(work[j] == 0.0);
      xfree(work);
      return;
}

SPM *spm_mul_mat(const SPM *A, const SPM *B)
{     /* multiply two sparse matrices (driver routine) */
      SPM *C;
      C = spm_mul_sym(A, B);
      spm_mul_num(C, A, B);
      return C;
}

PER *spm_create_per(int n)
{     /* create permutation matrix */
      PER *P;
      int k;
      xassert(n >= 0);
      P = xmalloc(sizeof(PER));
      P->n = n;
      P->row = xcalloc(1+n, sizeof(int));
      P->col = xcalloc(1+n, sizeof(int));
      /* initially it is identity matrix */
      for (k = 1; k <= n; k++)
         P->row[k] = P->col[k] = k;
      return P;
}

void spm_check_per(PER *P)
{     /* check permutation matrix for correctness */
      int i, j;
      xassert(P->n >= 0);
      for (i = 1; i <= P->n; i++)
      {  j = P->row[i];
         xassert(1 <= j && j <= P->n);
         xassert(P->col[j] == i);
      }
      return;
}

void spm_delete_per(PER *P)
{     /* delete permutation matrix */
      xfree(P->row);
      xfree(P->col);
      xfree(P);
      return;
}

/* eof */