Re: Non square matrix transpose?

["E. Robert Tisdale" <edwin at netwood dot net>]

Brian Gough wrote:

> E. Robert Tisdale writes:
> > Are you sure that there is no other missing functionality?  Or do
> > you mean you to just keep tacking functions onto the library as
> > people ask for them?
>
> The library is under development (version 0.x),
> so more functions will be added/removed/redesigned in the future.
>
> When we get to version 1.0 the API will be frozen
> (one of the criteria for the 1.0 release is that we are happy with the API).
> Any changes in the 1.x series will be backwards compatible.
>
> Of course, we might also start a 2.x version based on lessons learned
> from 1.0, and that would not have to be backwards compatible.
>
> Like most free software projects
> there is no a priori specification of the complete system --
> just a goal to work towards
> and some guidelines/conventions for getting there.

No.
A library is NOT just another free software project.
The API is a contract between you, the library developer,
and application programmers.
It specifies what (ANSI C) application programmers should write
to access the functionality of the GSL.
You shouldn't be feeling your way like a blind man with a white cane.
It isn't necessary.

You should just go ahead and specify all of the functionality
that application programmers might require.
Basically, the GSL should include vector and matrix functions
analogous to each of the scalar operators and functions
supported by the (ANSI C) programming language.
Your API fails if application programmers are obliged
to implement any missing functionality as a loop
over subscripted vectors and/or matrices.
These loops introduce too many opportunities to make mistakes
that are hard to detect and debug.
They are also inherently serial and preclude the possibility
of parallel implementation at some time in the future.
Don't worry about supporting all of the functionality
specified in your API.
All the GSL really needs to support is the functionality
required to implement the functionality specified by your API
so that application programmers can implement it themselves
and share their implementations as they see fit.

I'm glad that you have acknowledged the possibility
that your implementation of the GSL may change.
But I'm not sure that you realize yet,
that you have painted yourself into a corner in a way
that makes it very hard to change the implementation
without changing the API.  If you change the API,
all of the existing application programs that use the GSL
will need to be recoded or discarded.
At the very least, what you need to do right now
is to hide the actual data representation
of block, vector and matrix objects
and provide application programmers with GSL functions
to access private data members indirectly.
At least, that way, you can change the data representation
and the GSL functions which access them and redistribute the GSL
without breaking any existing application programs.

> >  It strikes me as just a bit inefficient to make a copy of a matrix
> > just so that you can reference the elements in a different order.
> > And, as I recall, the BLAS library will operate just as happily on
> > the transpose view of a matrix as upon the matrix itself.
>
> For BLAS operations one should certainly use the transpose argument
> and not physically transpose the matrix.  In other applications
> it may be necessary to rearrange the elements though.
>
> >  Why couldn't you just add a field of type CBLAS_TRANSPOSE_t to
> > your matrix struct which tells GSL library functions whether the
> > data array represents a matrix or it's transpose?
>
> Here we follow the existing NAG and IMSL matrix conventions,
> which are based on a fixed ordering
> and don't generally use a transpose argument
> (except in their own BLAS/LAPACK-type routines).

The existing NAG and IMSL  matrix conventions
are constrained by deficiencies in the Fortran 77 programming language.
Specifically, Fortran 77 provides no mechanism for allocating
temporary workspace dynamically.
It isn't practical for Fortran 77 subroutines to repack matrix objects
into representations which are more suitable to the algorithm
so the Fortran 77 application programmer is obliged to repack
the matrix before calling the Fortran 77 subroutine.
Again, this introduces opportunities for the application programmer
to make mistakes that are difficult to detect and debug.
The ANSI C programming language, on the other hand,
provides a convenient mechanism for allocating workspace
for a temporary matrix representation suitable for NAG and IMSL.
The GSL function should automatically repack the matrix,
if necessary, before calling the NAG or IMSL library routine.
If you like, the GSL function could issue a message (called a tip)
to the application programmer complaining about the need
to repack the matrix but the tip should probably be ignored
in most cases because application programmers can seldom
reorganize their code to avoid the repacking and, in general,
it is safer to let the GSL function repack it for them.

Please reconsider the Abstract Data Type (ADT)
that the GSL is supposed to support.
Arbitrary strides (negative and zero as well as positive)
occur naturally in numerical application programs.
Ignoring them only compels application programmers
to write loops to convert them to more restrictive representations.
This wastes storage space and increases the probability
the the application programmer will make a mistake.
The overall result is a decrease in performance and reliability.