Glas List: Re: [glas] concepts, algorithms etc

From: Matthias Troyer (troyer_at_[hidden])
Date: 2005-03-11 10:18:29

On Mar 11, 2005, at 6:16 AM, Karl Meerbergen wrote:
> My question is: why do we need a VectorSpace concept?

When writing the iterative eigensolver template library (IETL), a 
generalization of the ITL to eigenvalue problems, we encountered the 
fact that the concepts used in the ITL were not general enough and had 
to introduce a vector space comcepts.(see 
http://www.comp-phys.org/software/ietl/ to obtain the IETL library)

The problems that made us introduce a vector space concept were:

1.) creation of new vectors
2.) obtain the dimension of a vector space
3.) projecting vectors from a "representation" vector space to the real 
"physical" vector space.

Before explaining the above points let me quickly show what functions 
and types we introduced:

For a model VS of the concept vector space we defined the following 
types

vectorspace_traits<VS>::vector_type
vectorspace_traits<VS>::scalar_type
vectorspace_traits<VS>::size_type

and the following functions:

   typename ietl::vectorspace_traits<VS>::vector_type new_vector(const 
VS& vs)
   typename ietl::vectorspace_traits<VS>::size_type vec_dimension(const 
VS& vs)
   void project(typename ietl::vectorspace_traits<VS>::vector_type& v, 
const VS& vs)

Now let me explain why we needed these functions and concepts by 
comparing to what the ITL did:

1.) new_vector:

Here the ITL assumed that a new vector of type V could be created by 
calling a constructor taking the size of the vector, e.g. given the 
dimension d of the vector space it assumed that a new vector can be 
constructed by

V x(s);

In our applications this caused several problems, where I will just 
mention two typical problems:
a) to solve the eigenvalue problem for a 3D partial differential 
equation on a regular mesh we most naturally use a three-dimensional 
Blitz array (or boost::multi_array) as our data type. These 3D arrays 
form a vector space, but have no constructor of the above type.
b) for parallel vectors we sometimes want to pass information on how 
the vector should be distributed over the nodes in the constructor.

The new_vector(vs) function solves these problems.

2.) vec_dimension:

allows us to obtain the dimension of the vector space, without creating 
a vector.

3.) project:

came up when we realized that in some applications the actual vector 
space we are interested is actually a subspace of the vector space 
spanned by the data type, defined by some constraints. One example is 
again solving a partial differential equation, with fixed boundary 
conditions. The variables might again be stored in a 3D array, but 
where the values in the boundary layer are fixed. Solving such an 
eigenvalue problem by an iterative eigensolver requires starting from a 
random starting vector. Obviously we cannot just fill the entire vector 
with random numbers, since this would give us wrong values on the 
boundaries. The project() function can then be used to make sure that 
the boundary values are correct, thus "projecting" the vector from the 
vector space spanned by the data type (with arbitrary boundary values) 
to the desired vector space (with fixed boundary values).

A student of mine will soon start a term project to add these vector 
space features to the ITL library, so that we can also use the ITL on 
our data types. We'll keep you updated on our progress.

Matthias