svector.txt

/**************************************************************************\

MODULE: svector

SUMMARY:

Template class for dynamic-sized sparse vectors.

The declaration

   SVec<T> v;

creates a zero-length vector.  To grow this vector to length n, execute

   v.SetLength(n)

Note that this call does not allocate any space for the elements.  Instead,
elements are allocated as they are set.  Since this vector is intended to be
used for sparse vectors, the number of elements allocated in the vector should
be much less than the length.  As elements get allocated, the default
constructor for T is called to initialize the elements.

The current length of a vector is available as v.length().

The i-th vector element (counting from 0) is accessed as v[i].  If the
macro NTL_RANGE_CHECK is defined, code is emitted to test if 0 <= i <
v.length().  This check is not performed by default.  

For old-time FORTRAN programmers, the i-th vector element (counting
from 1) is accessed as v(i).

If v is const and a "zero" element is accessed (an element for which storage
has not been allocated), then a reference to a static const "zero" element is
returned.  To provide this reference, one must also declare somewhere

  template<> const T& zero_ref<T>() { return ...; }

If v is not const and an element that has not been allocated is accessed, 
space is allocated for the element and initialized to zero using clear().
Because of this automatic allocation of elements, care must be taken when
iterating through the elements of the vector.  The most efficient method
for iterating through the elements is to make use of the methods nvalues(),
indices(), and values().  Alternatively, cast v to const before attempting
to read all of the elements of the vector.

Let n = v.length().  Calling v.SetLength(m) with m <= n sets the
current length of v to m but does not call any destructors or free
any space.  As discussed above, calling v.SetLength(m) with m > n does not
allocate any new space or change any existing elements.  Space is allocated
strictly as needed.

v.MaxLength() is the largest value of n for which v.SetLength(n) was invoked.
v.SetMaxLength(n) only changes MaxLength if n is greater than the previous
value of MaxLength.  MaxLength is not of much use in svector.  It is provided
simply for compatibility with (non-sparse) vectors.

When v's destructor is called, all constructed elements will be destructed, 
and all space will be relinquished.

Space for storing elements is allocated by the underlying (non-sparse) vector
class.  See the vector module for details regarding space management.

Note that when new space is required to store the non-zero elements, the space
is reallocated using realloc, and thus the addresses of vector elements may
change, possibly creating dangling references to vector elements.  This
behaviour is also an issue with the (non-sparse) vector classes; however,
since new space may be allocated when elements of the sparse vector are simply
referenced (as opposed to when SetLength() is called), extra care must be
taken when passing references to vector elements around.

v.allocated() is the number of elements which have been allocated, which
indicates an upper bound on the number of non-zero elements in the vector.
The method v.SetAlloc(n) may be used to force allocation for up to n elements
in advance of setting values.

If the type T supports standard arithmetic operations such as add(x,a,b),
sub(x,a,b), mul(x,a,b), negate(x,a), and the operators +, -, and * (as most
NTL types do) then many standard math function can be declared using
NTL_math_svector_decl(T,vec_T,svec_T) and implemented with
NTL_math_svector_impl(T,vec_T,svec_T).



\**************************************************************************/

template <class T>
class SVec {
public:

  /* Default constructor.  Initial length is zero and no memory
   * is allocated.  */
  SVec(); 

  /* Copy constructor.  Allocates exactly the right amount of memory
   * to hold the contents of other and uses T's assignment operator to
   * do the copying.  If compact then zero values in other are not copied.
   * The "fixed" status of other is not copied.  */
  SVec(const SVec& other, bool compact = false);

  /* Assignment.  Performs an element-wise assignment using T's assignment
   * operator.  If this is "fixed", other must have the same length as 
   * this.  */
  SVec& operator=(const SVec& other);

  /* Initialize with a specific length and allocation.  If alloc is zero then
   * no memory is allocated until needed (and then the default allocation is
   * made).   */
  SVec(INIT_SIZE_TYPE, long length, long alloc = 0);

  /* Destructor.  Free all allocated memory.  */
  ~SVec();

  /* Release all allocated space and set to length 0.  */
  void kill();

  /* Ensure enough memory is allocated for n non-zero entries.  */
  void SetAlloc(long n);

  /* Set maximum length to n.  Does not allocate memory and does not change
   * the current length.  Here for compatibility with non-sparse vector
   * classes.  */
  void SetMaxLength(long n);

  /* Sets length to n and prohibits all future length changes.
   * FixLength may only be invoked immediately after the default
   * construction or kill.
   *
   * The kill operation is also subsequently prohibited, and swap is
   * allowed on fixed length vectors of the same length.
   *
   * FixLength is provided mainly to implement smat_T, to enforce
   * the restriction that all rows have the same length.  */
  void FixLength(long n);

  /* Set current length to n.  If the length is being set smaller than it was,
   * values at the end of the vector may disappear.  If n is being set larger,
   * existing values are not changed.   */
  void SetLength(long n);

  /* Current length.  */
  long length() const;

  /* Has this vector been fixed?  */
  bool fixed() const;

  /* Maximum length ever achieved.  Here for compatibility with non-sparse
   * vector classes.  */
  long MaxLength() const;

  /* The number of objects for which space has been allocated but not
   * necessarily initialized.  */
  long allocated() const;
    
  /* Indexing operation, starting from 0.  This version returns a non-const
   * reference to a T and will allocate memory for the vector element as
   * needed.  If you are just reading vector elements and you want to make
   * sure you maintain the sparsity of the vector, use the const version
   * instead.   */
  T& operator[](long i);

  /* Indexing operation, starting from 0.  This version can be used with a
   * const object and returns a const reference to a T.  No memory will be
   * allocated.  If the i'th vector element is not allocated, a const
   * reference to zero will be returned.   */
  const T& operator[](long i) const;

  /* Indexing operation, starting from 1.  See operator[] for details.  */
  T& operator()(long i);

  /* Indexing operation, starting from 1.  See operator[] for details.  */
  const T& operator()(long i) const;

  /* Direct access to arrays, use with care.  */
  long nvalues() const;
  const long* indices() const;
  T* values();
  const T* values() const;

  /* Indexing with no range checking.  See operator[] for details.  */
  T& RawGet(long i);
  const T& RawGet(long i) const;

  /* Returns position of a in the vector, or -1 if it is not there.
   * The search is conducted from position 0 to MaxAlloc()-1 of the vector,
   * and an error is raised if the object is found at position MaxLength()
   * or higher (in which case a references an uninitialized object).
   * Note that if NTL_CLEAN_PTR flag is set, this routine takes
   * linear time, and otherwise, it takes constant time.  */
  long position(const T& a) const;

  /* Returns position of a in the vector, or -1 if it is not there.
   * The search is conducted from position 0 to length()-1 of the vector.
   * Note that if NTL_CLEAN_PTR flag is set, this routine takes
   * linear time, and otherwise, it takes constant time.  */
  long position1(const T& a) const;

  /* Eliminate zero values from vector to speed up read-only access.  */
  void compact();

  /* Set all elements to zero.  Does not change length or release any 
   * memory.  */
  void clear();

  /* Swaps this and other by swapping internal pointers.  If either this or
   * other is fixed, then both must be have the same length.  The "fixed"
   * status of the vectors is not swapped.  */
  void swap(SVec& other);
};



/**************************************************************************\

                       Some utility routines

\**************************************************************************/

// test if x is the zero vector
long IsZero(const SVec<T>& x);

// make x the zero vector (without changing its length)
void clear(SVec<T>& x);

// swap x and y by swapping pointers
void swap(SVec<T>& x, SVec<T>& y);

// convert between Vec<T> and SVec<T>
void conv(SVec<T>& dest, const Vec<T>& src);
void conv(Vec<T>& dest, const SVec<T>& src);

// create copy of a with length exactly n (input is truncated or padded
// with zeros as necessary)
void VectorCopy(SVec<T>& x, const SVec<T>& a, long n);

// create copy of a with length exactly n (input is truncated or padded
// with zeros as necessary)
SVec<T> VectorCopy(const SVec<T>& a, long n);



/**************************************************************************\

                              Equality Testing

The tests are performed using the underlying operator == for T.

\**************************************************************************/

long operator==(const SVec<T>& a, const SVec<T>& b);
long operator!=(const SVec<T>& a, const SVec<T>& b);



/**************************************************************************\

                             Input/Output

Elements are written and read using the underlying I/O operators << and >>
for T.

The I/O format for a sparse vector v with length n is:

   <i0 v0 i1 v1 ... n>

where the i's are indices and the v's are (non-zero) values.  The i's are in
order and 0 <= i < n.

The method OutputVector() will write out the sparse vector as a dense vector
compatible with Vec<T>.

   [v[0] v[1] ... v[n-1]]

This will waste storage if the vector is really sparse, but can be read back
as a dense vector.  

The input method will read either format.

\**************************************************************************/

istream& operator>>(istream&, SVec<T>&);
ostream& operator<<(ostream&, const SVec<T>&);
ostream& OutputVector(ostream&, const SVec<T>&);



/************************************************************************** \

                             Arithmetic

The arithmetic functions and operators can be declared with 
NTL_math_svector_decl(T,vec_T,svec_T) and implemented with 
NTL_math_svector_impl(T,vec_T,svec_T).

These methods require add(x,a,b), sub(x,a,b), mul(x,a,b), negate(x,a),
and the +, -, * operators for the underlying type T.

\**************************************************************************/

void mul(svec_T& x, const svec_T& a, const T& b);
void mul(svec_T& x, const T& a, const svec_T& b);

void add(svec_T& x, const svec_T& a, const svec_T& b);
void sub(svec_T& x, const svec_T& a, const svec_T& b);
void negate(svec_T& x, const svec_T& a);

void InnerProduct(T& x, const svec_T& a, const svec_T& b);

svec_T operator+(const svec_T& a, const svec_T& b);
svec_T operator-(const svec_T& a, const svec_T& b);
svec_T operator-(const svec_T& a);

svec_T operator*(const svec_T& a, const T& b);
svec_T operator*(const T& a, const svec_T& b);

T operator*(const svec_T& a, const svec_T& b);

svec_T& operator+=(svec_T& x, const svec_T& a);
svec_T& operator-=(svec_T& x, const svec_T& a);


/* mixed vector/svector math methods */

void mul(vec_T& x, const svec_T& a, const T& b);
void mul(vec_T& x, const T& a, const svec_T& b);

void add(vec_T& x, const vec_T& a, const svec_T& b);
void add(vec_T& x, const svec_T& a, const vec_T& b);
void sub(vec_T& x, const vec_T& a, const svec_T& b);
void sub(vec_T& x, const svec_T& a, const vec_T& b);
void negate(vec_T& x, const svec_T& a);

void InnerProduct(T& x, const vec_T& a, const svec_T& b);
void InnerProduct(T& x, const svec_T& a, const vec_T& b);

vec_T operator+(const vec_T& a, const svec_T& b);
vec_T operator+(const svec_T& a, const vec_T& b);
vec_T operator-(const vec_T& a, const svec_T& b);
vec_T operator-(const svec_T& a, const vec_T& b);

T operator*(const vec_T& a, const svec_T& b);
T operator*(const svec_T& a, const vec_T& b);

vec_T& operator+=(vec_T& x, const svec_T& a);
vec_T& operator-=(vec_T& x, const svec_T& a);