alpar@1: /* glpios04.c (operations on sparse vectors) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpios.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_create_vec - create sparse vector
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  IOSVEC *ios_create_vec(int n);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_create_vec creates a sparse vector of dimension n,
alpar@1: *  which initially is a null vector.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine returns a pointer to the vector created. */
alpar@1: 
alpar@1: IOSVEC *ios_create_vec(int n)
alpar@1: {     IOSVEC *v;
alpar@1:       xassert(n >= 0);
alpar@1:       v = xmalloc(sizeof(IOSVEC));
alpar@1:       v->n = n;
alpar@1:       v->nnz = 0;
alpar@1:       v->pos = xcalloc(1+n, sizeof(int));
alpar@1:       memset(&v->pos[1], 0, n * sizeof(int));
alpar@1:       v->ind = xcalloc(1+n, sizeof(int));
alpar@1:       v->val = xcalloc(1+n, sizeof(double));
alpar@1:       return v;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_check_vec - check that sparse vector has correct representation
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_check_vec(IOSVEC *v);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_check_vec checks that a sparse vector specified by
alpar@1: *  the parameter v has correct representation.
alpar@1: *
alpar@1: *  NOTE
alpar@1: *
alpar@1: *  Complexity of this operation is O(n). */
alpar@1: 
alpar@1: void ios_check_vec(IOSVEC *v)
alpar@1: {     int j, k, nnz;
alpar@1:       xassert(v->n >= 0);
alpar@1:       nnz = 0;
alpar@1:       for (j = v->n; j >= 1; j--)
alpar@1:       {  k = v->pos[j];
alpar@1:          xassert(0 <= k && k <= v->nnz);
alpar@1:          if (k != 0)
alpar@1:          {  xassert(v->ind[k] == j);
alpar@1:             nnz++;
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(v->nnz == nnz);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_get_vj - retrieve component of sparse vector
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  double ios_get_vj(IOSVEC *v, int j);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine ios_get_vj returns j-th component of a sparse vector
alpar@1: *  specified by the parameter v. */
alpar@1: 
alpar@1: double ios_get_vj(IOSVEC *v, int j)
alpar@1: {     int k;
alpar@1:       xassert(1 <= j && j <= v->n);
alpar@1:       k = v->pos[j];
alpar@1:       xassert(0 <= k && k <= v->nnz);
alpar@1:       return (k == 0 ? 0.0 : v->val[k]);
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_set_vj - set/change component of sparse vector
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_set_vj(IOSVEC *v, int j, double val);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_set_vj assigns val to j-th component of a sparse
alpar@1: *  vector specified by the parameter v. */
alpar@1: 
alpar@1: void ios_set_vj(IOSVEC *v, int j, double val)
alpar@1: {     int k;
alpar@1:       xassert(1 <= j && j <= v->n);
alpar@1:       k = v->pos[j];
alpar@1:       if (val == 0.0)
alpar@1:       {  if (k != 0)
alpar@1:          {  /* remove j-th component */
alpar@1:             v->pos[j] = 0;
alpar@1:             if (k < v->nnz)
alpar@1:             {  v->pos[v->ind[v->nnz]] = k;
alpar@1:                v->ind[k] = v->ind[v->nnz];
alpar@1:                v->val[k] = v->val[v->nnz];
alpar@1:             }
alpar@1:             v->nnz--;
alpar@1:          }
alpar@1:       }
alpar@1:       else
alpar@1:       {  if (k == 0)
alpar@1:          {  /* create j-th component */
alpar@1:             k = ++(v->nnz);
alpar@1:             v->pos[j] = k;
alpar@1:             v->ind[k] = j;
alpar@1:          }
alpar@1:          v->val[k] = val;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_clear_vec - set all components of sparse vector to zero
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_clear_vec(IOSVEC *v);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_clear_vec sets all components of a sparse vector
alpar@1: *  specified by the parameter v to zero. */
alpar@1: 
alpar@1: void ios_clear_vec(IOSVEC *v)
alpar@1: {     int k;
alpar@1:       for (k = 1; k <= v->nnz; k++)
alpar@1:          v->pos[v->ind[k]] = 0;
alpar@1:       v->nnz = 0;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_clean_vec - remove zero or small components from sparse vector
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_clean_vec(IOSVEC *v, double eps);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_clean_vec removes zero components and components
alpar@1: *  whose magnitude is less than eps from a sparse vector specified by
alpar@1: *  the parameter v. If eps is 0.0, only zero components are removed. */
alpar@1: 
alpar@1: void ios_clean_vec(IOSVEC *v, double eps)
alpar@1: {     int k, nnz;
alpar@1:       nnz = 0;
alpar@1:       for (k = 1; k <= v->nnz; k++)
alpar@1:       {  if (fabs(v->val[k]) == 0.0 || fabs(v->val[k]) < eps)
alpar@1:          {  /* remove component */
alpar@1:             v->pos[v->ind[k]] = 0;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* keep component */
alpar@1:             nnz++;
alpar@1:             v->pos[v->ind[k]] = nnz;
alpar@1:             v->ind[nnz] = v->ind[k];
alpar@1:             v->val[nnz] = v->val[k];
alpar@1:          }
alpar@1:       }
alpar@1:       v->nnz = nnz;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_copy_vec - copy sparse vector (x := y)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_copy_vec(IOSVEC *x, IOSVEC *y);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_copy_vec copies a sparse vector specified by the
alpar@1: *  parameter y to a sparse vector specified by the parameter x. */
alpar@1: 
alpar@1: void ios_copy_vec(IOSVEC *x, IOSVEC *y)
alpar@1: {     int j;
alpar@1:       xassert(x != y);
alpar@1:       xassert(x->n == y->n);
alpar@1:       ios_clear_vec(x);
alpar@1:       x->nnz = y->nnz;
alpar@1:       memcpy(&x->ind[1], &y->ind[1], x->nnz * sizeof(int));
alpar@1:       memcpy(&x->val[1], &y->val[1], x->nnz * sizeof(double));
alpar@1:       for (j = 1; j <= x->nnz; j++)
alpar@1:          x->pos[x->ind[j]] = j;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_linear_comb - compute linear combination (x := x + a * y)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_linear_comb(IOSVEC *x, double a, IOSVEC *y);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_linear_comb computes the linear combination
alpar@1: *
alpar@1: *     x := x + a * y,
alpar@1: *
alpar@1: *  where x and y are sparse vectors, a is a scalar. */
alpar@1: 
alpar@1: void ios_linear_comb(IOSVEC *x, double a, IOSVEC *y)
alpar@1: {     int j, k;
alpar@1:       double xj, yj;
alpar@1:       xassert(x != y);
alpar@1:       xassert(x->n == y->n);
alpar@1:       for (k = 1; k <= y->nnz; k++)
alpar@1:       {  j = y->ind[k];
alpar@1:          xj = ios_get_vj(x, j);
alpar@1:          yj = y->val[k];
alpar@1:          ios_set_vj(x, j, xj + a * yj);
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  ios_delete_vec - delete sparse vector
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpios.h"
alpar@1: *  void ios_delete_vec(IOSVEC *v);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine ios_delete_vec deletes a sparse vector specified by the
alpar@1: *  parameter v freeing all the memory allocated to this object. */
alpar@1: 
alpar@1: void ios_delete_vec(IOSVEC *v)
alpar@1: {     /* delete sparse vector */
alpar@1:       xfree(v->pos);
alpar@1:       xfree(v->ind);
alpar@1:       xfree(v->val);
alpar@1:       xfree(v);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /* eof */