First let me clarify how valarray<> is used in may C++ dense linear
algebra class library. When I was planning the development of LinAlgPack
I wanted to use and implementation that would allow me to directly call
BLAS operations because I thought that that was my best shot at getting
the best performance. From what I know of valarray<> and slice<> these
fit the bill perfectly and I just assumed that the elements in valarray<>
would be stored sequentially. With 32 bit operating systems and flat
address spaces why would you do it any other way. However, for Level-1
BLAS operations I thought that there was a chance that optimized
implementations of valarray<> would outperform BLAS on these platforms.
When I designed VectorSlice I had to decide on whether it should provide a
view of a raw C array or a valarray<> object. By basing VectorSlice view
on a C array you could open up the possibility of allowing clients to
allocate the arrays themselfs, possibly in other languages like Fortran,
and then through a C interface they would have access to your entire
linear algebra library and any other specialized algorithms you have
written yourself in terms of VectorSlice. By using valarray<> as the
basis however you would have to copy all data to and from valarray<>
perform the same tasks. On the other hand by using valarray<> you could
validate that you never create a view into data that does not exist and
you have the possibility of using valarray<> to implement any of your
linear algebra operations if it is adventagous.
For example, if you create a VectorSlice using the constructor:
VectorSlice(double* p, size_t size, ptr_diff_t stride = 1)
how do you know if:
(1) p even points to an allocated array
(2) p + size * abs(stride) is not past the end of the allocated array
However if you basis the view on a valarray<> (or any other container for
that matter then with the constructor:
VectorSlice(valarray<double>& v, size_t start, size_t size, ptr_diff_t
stride = 1)
then you you can validate (1) and (2) above very easily.
I used valarray<> for the basis of my VectorSlice class because of the
built in validation shown above and because if core of your application is
in C++ and calls other legacy code then all (or at least most) data
allocation will occur in
C++ and the disadvantages of using valarray<> are greatly diminished.
In reallity the standard library user never even sees or works with a
valarray. Heres why.
Here is a simplified declaration for VectorSlice:
class VectorSlice {
public:
// public types
typedef ... value_type;
typedef ... size_type;
typedef ... difference_type;
typedef stride_iter<value_type*
,difference_type, ... > iterator;
typedef ... const_iterator;
typedef valarray<value_type> valarray;
// constructors
VectorSlice(valarray& v, size_type start, size_type size
, differenc_type stride = 1);
// size
size_type size() const;
// assignment (copies to target valarray)
VectorSlice& operator=(const VectorSlice& vs_rhs);
VectorSlice& operator=(value_type alpha);
// iterator access
iterator begin();
iterator end();
const_iterator begin() const;
const_iterator end() const;
// subregion access
VectorSlice operator()(size_type lbound, size_type ubound);
VectorSlice operator()(Range1D rng);
const VectorSlice operator()(size_type lbound, size_type ubound) const;
const VectorSlice operator()(Range1D rng) const;
// representation access (this will never change)
value_type* start_ptr();
const value_type* start_ptr() const;
difference_type stride() const;
// implementation access (only for implementators
// , breach of encapsulation)
valarray& v();
const valarray& v() const;
}; // end class VectorSlice
The library user never creates valarray<> objects directly, instead they
create Vector objects:
class Vector {
// public types (the same as for VectorSlice)
...
// constructors (all of the standard valarray<> constructors
Vector(size_type size);
...
// VectorSlice public interface
...
// Conversion to a VectorSlice
operator VectorSlice();
operator const VectorSlice() const;
private:
valarray v_; // no storage over head over a valarray<>
}; // end class Vector
Then when a user needs a VectorSlice to call a function they create a
Vector object then perform a conversion to a VectorSlice object. Here is
an example:
// Of course you would use BLAS xDOT to do this in a real implementation
inline value_type norm_2(const VectorSlice& vs) {
value_type result = 0.0;
for(VectorSlice::const_iterator it =vs.begin(); itr!=vs.end();++itr)
result += (*itr) * (*itr);
return sqrt(result);
}
void main() {
Vector v(10);
v = 1.0;
cout << dot(v);
}
In the above example the user never even know that a valarray<> was being
used and the implementation for dot() did not either. If fact nothing was
exposed about the implementation. Hell, Vector could have been a sparse
vector implemented as a linked list and the iterator type could have been
as STL linked list iterator.
My GenMatrix and GenMatrixSlice classes have the same roles as Vector and
VectorSlice.
>
> The library developer SHOULD hide the actual data representation
> and provide implementation independent methods
> which the portable application programmer can use to access data.
> There is no way to deny the application program direct access
> to the underlying data representation but the library developer
> is not obliged to support application programs which do so.
>
Does anyone think that it is unreasonable to say that:
"The data elements that VectorSlice provides a view into are, and forever
shall be stored in a continous array separated by stride possitions."
Fortran and the BLAS have gotten by with this for decades. This brings up
another point. Why do we even bother writting C++ class libraries for
dense linear algebra when the data stuctures we are abstracting are the
simplest you can imagine? Is it really encapsulation we are going for or
is it abstraction?
For example, compare these two matrix-vector multiplication functions for:
y = alpha * A * x + beta * y
// Raw data implemtation
void mult_1(double alpha, const double A[], size_t m, size_t n, size_t lda
, const double x[], double beta, double y[])
{
// use loops, blocking or call the BLAS, whatever you want.
...
}
// Matrix / Vector abstraction implementation
void mult_2(double alpha, const GenMatrixSlice& A, const VectorSlice& x
, double beta, VectorSlice* y)
{
// use loops with iterators or raw pointers, valarray<>
// (not recommended for the general user), or call the BLAS.
...
}
Why would we prefere mult_2() to mult_1() if we have the chooce?
(1) There are lesss arguments to keep track of and are grouped in a
logical way.
(2) We can build in validation within mult_2() to ensure that A, x and y
are of the correct sizes for debugging at no extra inconviance to the
user.
Neither (1) or (2) has anything to do with encapsulation, only
abstraction. Given that I have a GenMatrixSlice A and VectorSlice objects
x and y I should be able to call functions of the form of mult_1() with
some assurance that this will be supported and portable. My god! Is
really unreasonable? I could not use a C++ class library if this were not
the case.
> It is important to distinguish the role played by the library developer
> from the role played by the application programmer.
> So far, Roscoe has been playing both roles and he wants to continue doing so.
> That's fine as long as he is the only one using his library.
> If he decided to change the underlying representation,
> he could just change all of his applications as well.
> But I hope, if he has convinced someone else to use his library,
> he would consider what effect changes in his library would have
> on their applications and devise some scheme to ensure
> that the impact would be minimized.
As I stated above, the general user never sees valarray<> and only creates
VectorSlice objects from Vector objects. I could remove valarray<> and
not have to change one line of code. I would only have to recompile it.
>
> The actual data representation is irrelevant to the application programmers
> if the library is complete and well designed but if they discover
> a deficiency in the library, they may require direct access to fix it.
> In this case, the application programmer switches hats
> and becomes a "library co-developer"
> in an implicit partnership with the original library developer.
> Usually, application programmers can get all the cooperation they need
> from the library developer as long as they don't require any additional
> commitments from the library developer.
>
The library user should be able to neglect the libraries implementation of
an operation and implement his or her own if they feel like it and be
assured the it will still be applicable in the next version of the
library. Remember, we are talkine about dense linear algebra here where
the most complicated data structure is:
double A[], size_t m, size_t n, size_t lda
If anyone is interested I have set of standards that I have been working
on for my own libraries. I have put a lot of work and thought into them
and will distribute the html files for them if desired.
(1) Naming sceme for linear algebra operations supported by functions that
is: a) consistant, b) compact, c) easy to translate from the mathematical
notation to the function identifier name.
(2) Template class specifications (compile time polymorphism) for a)
sparse unstored vectors of the form (val,i) and b) sparse unstorted
coordinate matrices of the form (val,i,j).
(3) Abstract interface to matices and Level-2 and 3 BLAS operations using
virtual functions to used when runtime polymorphism is needed.
Matrix (rows, cols)
|
MatrixWithOp (y = alpha*op(A)*x + beta*y, Y = alpha*op(A)*op(B) + beta*Y)
|
MatrixWithOpFactorized (y = alpha * inv(op(A)) * x, ... )
--------------------------------------
| Roscoe Bartlett |
| Ph.D. student |
| Department of Chemical Engineering |
| Carnegie Mellon University |
| email: roscoe@andrew.cmu.edu |
--------------------------------------