ode-initval implicit solvers and development

ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

my work with GSL ODE solvers is progressing (slowly, I am doing this on 
my own free time). I've written a modified Newton iteration method to 
solve the non-linear equation system that come up with implicit 
Runge-Kutta methods. Three very basic implicit solvers included use the 
solver. Also, I've made some changes and additions to ode-initval tests. 
Those are the highlights. The current work on gsl-1.11/ode-initval 
directory can be found here:

http://iki.fi/tuomo.keskitalo/temp/gsl-1.11-ode-initval-additions-1.tar.gz

Please note the comments on file Changes.txt. Any comments you might 
have are welcome!

Fifth order implicit Radau method would be nice to add maybe later on 
(even these new implicit solvers do not seem to be able to cope with 
highly stiff, badly scaled problems), but I think that I will continue 
next with a multistep framework and some of the methods in e.g. the VODE 
package. More tests are also in the queue.

Regards,
Tuomo

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Brian Gough]

At Tue, 30 Sep 2008 20:06:49 +0300,
Tuomo Keskitalo wrote:
> 
> http://iki.fi/tuomo.keskitalo/temp/gsl-1.11-ode-initval-additions-1.tar.gz
> 
> Please note the comments on file Changes.txt. Any comments you might 
> have are welcome!

Great! I will start looking at it...

Re: ode-initval implicit solvers and development

[Brian Gough]

At Tue, 30 Sep 2008 20:06:49 +0300,
Tuomo Keskitalo wrote:
> my work with GSL ODE solvers is progressing (slowly, I am doing this on 
> my own free time). I've written a modified Newton iteration method to 
> solve the non-linear equation system that come up with implicit 
> Runge-Kutta methods. Three very basic implicit solvers included use the 
> solver. Also, I've made some changes and additions to ode-initval tests. 

This is good, I would like to replace the existing rkimp solvers with
these in a future release.  The main outstanding change needed is to
avoid recomputing the jacobian J(t,x) in the second call to
modnewton1_init.  The parameters t and x are the same so it computed
once in the caller and passed to modnewton1_init as a parameter.

-- 
Brian Gough

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

while waiting for a reference book on multistep methods to arrive, I've 
made some more modifications to current ode-initval, which I hope are 
now close to being ready to be included into GSL. The .git directory 
that should include the changes can be found at

http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git

The changes to current GSL git repository include:

- New implicit solvers: imprk4 (replaces imprk42 in previous test
package), imprk2 (replaces imprk21 in previous test package) and
impeuler, which use modified Newton iteration to solve the implicit
equations. I decided to dump imprk42 and imprk21 because I was unsure
how to calculate the error estimates for them. I think that these new
implicit methods could be used instead of rk4imp, rk2imp, gear2 and
gear1, which use functional iteration (which is inefficient for stiff
problems.) However, even these new solvers can not be used efficiently
for truly nasty stiff problems, and better solvers are still needed.

- modnewton1.c: a modified Newton solver for solution of non-linear
equations of implicit solvers. Jacobian evaluation has been moved from
modnewton1 to imprk2, imprk4 and impeuler.

- evolve.c: modification to decrease step size if stepper fails,
instead of giving up immediately.

- evolve.c: A bug fix: Exit with GSL_FAILURE if step size reaches
machine precision instead of continuing with old step size.

- test.c: added Robertson stiff test

- test.c: embedded test_oregonator to test_compare_stiff_probelms to
test several stiff problems. Note: the tested integration interval is
rather short due to inefficiency of solvers other than bsimp in these
problems.

- test.c: modified test_compare_vanderpol and
test_compare_stiff_problems to check against results from first
ode-solver. Changed rk4 or bsimp to first place.

- test.c: modified stepfn and stepfn2 and expanded these to be tested
with all explicit solvers

- test.c: removed stepfn3 and added test_broken and test_stepsize_fail
to check that evolve decreases step size below machine precision and
exits with a failure code in the end if user functions fail
continuously.

Changes in ode-initval.texi:

- Removed untrue sentence concerning stepping functions: "The
step-size @var{h} will be set to the step-size which caused the
error."

- Changed gsl_odeiv_system variable name from dydt to sys for clarity.

- Added "explicit" or "implicit" to stepper introductions for clarity.

- Reformulated description of evolve_apply and step_apply and made
changes to reflect the modifications to the code.

- Clarified the point on integrating over discontinuities. It is best
to integrate in sequences.

- Removed the text suggesting the user to force a step size to
integrate over the edge of a discontinuity. The applicability of this
kind of a numerical tweak depends on the case, and should in my
opinion not be included in a reference book.

- Added literature references.

Please feel free to comment on the proposed changes.

Regards,
Tuomo

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Brian Gough]

At Sun, 02 Nov 2008 19:34:53 +0200,
Tuomo Keskitalo wrote:
> while waiting for a reference book on multistep methods to arrive, I've 
> made some more modifications to current ode-initval, which I hope are 
> now close to being ready to be included into GSL. The .git directory 
> that should include the changes can be found at
> 
> http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git

Thanks, I've retrieved it successfully (much easier than patches by email).

For those new to git the relevant commands to download it are

git remote add tuomo http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git
git fetch tuomo 

The branch name is then tuomo/master.

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

I've been working on an implementation of the Adams linear multistep 
methods with functional iteration for GSL ode-initval (msadams), and 
finally got a first beta for you to test and comment on. The source is 
available at the git repository mentioned below.

Please note the comments on file ode-initval/ChangeLog-Keskitalo. The 
biggest change is an addition to the step and control object structures 
to allow passing error tolerance to steppers. I'm afraid a change like 
this is really needed.

I'm planning to implement also the backward differentiation formula 
(BDF) methods with Newton iteration for stiff problems.

Kind regards,
Tuomo

On 11/03/2008 08:08 PM, Brian Gough wrote:

> At Sun, 02 Nov 2008 19:34:53 +0200,
> Tuomo Keskitalo wrote:
>> while waiting for a reference book on multistep methods to arrive, I've 
>> made some more modifications to current ode-initval, which I hope are 
>> now close to being ready to be included into GSL. The .git directory 
>> that should include the changes can be found at
>>
>> http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git
> 
> Thanks, I've retrieved it successfully (much easier than patches by email).
> 
> For those new to git the relevant commands to download it are
> 
> git remote add tuomo http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git
> git fetch tuomo 
> 
> The branch name is then tuomo/master.

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Brian Gough]

At Sat, 24 Jan 2009 13:52:22 +0200,
Tuomo Keskitalo wrote:
> Please note the comments on file ode-initval/ChangeLog-Keskitalo. The 
> biggest change is an addition to the step and control object structures 
> to allow passing error tolerance to steppers. I'm afraid a change like 
> this is really needed.

I had a look at the code (which all looks pretty sensible) and have
a couple of questions.

- Currently the control object is only used by the multistep adams
stepper.  What other types of methods might need access to the error
control?

- In modnewton1.c there is the comment

    /* Stopping tolerance for Newton iteration */
    /* FIXME: tol should be user defined */

How would you handle the newton tolerance--do we need a way to pass
that in as well?

Valgrind detects an uninitialised variable in msadams.c, which is
flagged with -Wall as well.


-- 
Brian Gough

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

On 02/01/2009 05:42 PM, Brian Gough wrote:

> - Currently the control object is only used by the multistep adams
> stepper.  What other types of methods might need access to the error
> control?

All steppers that do some kind of iterative solution internally. For 
example the methods that use Newton iteration.

> - In modnewton1.c there is the comment
> 
>     /* Stopping tolerance for Newton iteration */
>     /* FIXME: tol should be user defined */
> 
> How would you handle the newton tolerance--do we need a way to pass
> that in as well?

Yes, I think that is needed. This would be done similarly to msadams, 
e.g. modnewton1 asks the tolerance via gsl_odeiv_control_errlev.

I plan to implement the Newton iteration for BDF methods first, and then 
modify modnewton1. I'd also like to implement the implicit 5th order 
Radau method, which also needs Newton iteration. This will take some 
time. It might be best to wait until then to see what other changes I'll 
have to make to the ode-initval interfaces. Hopefully none.

> Valgrind detects an uninitialised variable in msadams.c, which is
> flagged with -Wall as well.

I'll add initialization to those, to get rid of the warnings, in the 
next release. Thanks for your comments!

I'd appreciate if anyone has the chance to test and report how msadams 
compares to other steppers in practice and if there are any problems 
with it!

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

On 02/02/2009 07:05 PM, Tuomo Keskitalo wrote:

> I plan to implement the Newton iteration for BDF methods first, and then 
> modify modnewton1. I'd also like to implement the implicit 5th order 
> Radau method, which also needs Newton iteration. This will take some 
> time. It might be best to wait until then to see what other changes I'll 
> have to make to the ode-initval interfaces. Hopefully none.

Ok, I've made a new beta release which includes new msbdf stepper and 
modifications to modnewton1 and the related steppers, among other 
changes. Please look at ode-initval/ChangeLog-Keskitalo for more 
information. Comments are welcome.

git source can be cloned to current directory with command
git clone http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git

> I'd appreciate if anyone has the chance to test and report how msadams 
> compares to other steppers in practice and if there are any problems 
> with it!

This still applies. msadams and msbdf are somewhat complex steppers 
(multistep & multiorder, Nordsieck form) and need to be tested.

I think that if I later on have time to implement Radau II A, it will 
not introduce any further changes to the framework. Stability detection 
should still be added to msbdf, otherwise it is feature-complete. 
Thereby, I hope to hear your opinion about this development branch. Can 
it be incorporated into GSL, or should it be e.g. an extension library?

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Brian Gough]

At Sun, 01 Mar 2009 16:37:41 +0200,
Tuomo Keskitalo wrote:
> git source can be cloned to current directory with command
> git clone http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git

Hi,
Do you need to run git-update-index in that directory? I don't see any
new entries in the log:

(16:33)$ git log
commit 2291c2b7575055d050b904d04bf7fd557dfd6fc0
Merge: a0d8bf4... 1c51c98...
Author: Tuomo Keskitalo <tkeskita@pulmunen.(none)>
Date:   Sat Jan 24 11:52:02 2009 +0200

    Merge branch 'master' of git://git.savannah.gnu.org/gsl

-- 
Brian Gough

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

On 03/03/2009 06:34 PM, Brian Gough wrote:

> At Sun, 01 Mar 2009 16:37:41 +0200,
> Tuomo Keskitalo wrote:
>> git source can be cloned to current directory with command
>> git clone http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git
> 
> Hi,
> Do you need to run git-update-index in that directory? I don't see any
> new entries in the log:

Sorry, I actually forgot to run git-update-server-info. Should be fixed now.

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Heikki Orsila]

On Thu, Mar 05, 2009 at 09:47:24PM +0200, Tuomo Keskitalo wrote:
> On 03/03/2009 06:34 PM, Brian Gough wrote:
> 
> >At Sun, 01 Mar 2009 16:37:41 +0200,
> >Tuomo Keskitalo wrote:
> >>git source can be cloned to current directory with command
> >>git clone http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git
> >
> >Hi,
> >Do you need to run git-update-index in that directory? I don't see any
> >new entries in the log:
> 
> Sorry, I actually forgot to run git-update-server-info. Should be fixed now.

You can do "chmod a+rx hooks/post-update" to run "git 
update-server-info" automatically on "git push". man githooks.

-- 
Heikki Orsila
heikki.orsila@iki.fi
http://www.iki.fi/shd

Re: ode-initval implicit solvers and development

[Brian Gough]

At Thu, 05 Mar 2009 21:47:24 +0200,
Tuomo Keskitalo wrote:
> 
> Sorry, I actually forgot to run git-update-server-info. Should be fixed now.

I've downloaded it successfully now, that was the command I should
have said--I got them mixed up.  I will look at it over the next few
days.

-- 
Brian Gough

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

On 03/01/2009 04:37 PM, Tuomo Keskitalo wrote:

> I think that if I later on have time to implement Radau II A, it will 
> not introduce any further changes to the framework. Stability detection 
> should still be added to msbdf, otherwise it is feature-complete. 
> Thereby, I hope to hear your opinion about this development branch. Can 
> it be incorporated into GSL, or should it be e.g. an extension library?

Now a work version of a heuristic stability enhancement has been added 
to msbdf, among other finetuning and corrections.

I've now also added function benchmark_precision to test.c for 
benchmarking. Apparently msbdf works badly for rhs_func_exp, but 
currently I think that it is a feature and not a bug. More info at 
ChangeLog-Keskitalo.

Git source at http://iki.fi/tuomo.keskitalo/gsl/ode-initval-additions/.git

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Tuomo Keskitalo]

Hello,

On 03/01/2009 04:37 PM, Tuomo Keskitalo wrote:

> Thereby, I hope to hear your opinion about this development branch. Can 
> it be incorporated into GSL, or should it be e.g. an extension library?

I've just merged changes from master and updated ode-initval 
documentation etc. I think that my ode-initval branch is about complete, 
for now.

What do you think about this? How do we proceed? I don't think that I 
will add any new functionality to my branch until I hear about the 
continuation.

Best regards,
Tuomo

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: ode-initval implicit solvers and development

[Brian Gough]

At Fri, 01 May 2009 17:06:09 +0300,
Tuomo Keskitalo wrote:
> I've just merged changes from master and updated ode-initval 
> documentation etc. I think that my ode-initval branch is about complete, 
> for now.

Great.

> What do you think about this? How do we proceed? I don't think that I 
> will add any new functionality to my branch until I hear about the 
> continuation.

I'll study the new changes when I get back to the UK on Wednesday.

Re: ode-initval implicit solvers and development

[Brian Gough]

At Fri, 01 May 2009 17:06:09 +0300,
Tuomo Keskitalo wrote:
> What do you think about this? How do we proceed? I don't think that I 
> will add any new functionality to my branch until I hear about the 
> continuation.

As this involves a change to the API, namely in the gsl_odeiv_control
and gsl_odeiv_step_type struct it will need to be in a new major
version.

I'll make the 1.13 release next month, with bug fixes for 1.12 only,
and then incorporate the new odes and other changes so that we can go
to a 2.0 release.

Since the new implicit RK routines are more robust that the old ones,
I suggest we just update them in place so that we don't need the
alternate name imprk--just rkimp.

I put some minor changes in the branch origin/tuomo at
http://src.network-theory.com/gsl.git

-- 
Brian Gough

GSL 2.0 roadmap

[Tuomo Keskitalo]

On 05/08/2009 01:51 PM, Brian Gough wrote:

[was: Re: ode-initval implicit solvers and development]

> As this involves a change to the API, namely in the gsl_odeiv_control
> and gsl_odeiv_step_type struct it will need to be in a new major
> version.
> 
> I'll make the 1.13 release next month, with bug fixes for 1.12 only,
> and then incorporate the new odes and other changes so that we can go
> to a 2.0 release.

What is the (approximate) roadmap and milestones to 2.0?

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: GSL 2.0 roadmap

[Brian Gough]

Sorry for the delay in replying.  Here's a list of the changes from
the NOTES file and my own org file.  I plan to make a branch for
release 2.0 after the next release 1.13 which will hopefully be next
week.  The 1.x branch should probably be maintained in parallel with
any bug fixes for a release or two but hopefully those are just minor
things now which we can easily merge.

* Changes for Release 2.0
Break binary compatibility, but keep source compatibility.
** Add a 'void *' to all workspaces, to allow for future changes.
** Disable deprecated functions
** Fix up the workspace_alloc functions so they have consistent names
   (add functions where needed, don't remove)
** Standardize function names, in particular VERB vs NOUN (e.g. _invert
vs _inverse). Also adopt a convection for functions which can operate
in place vs use of workspace (e.g linalg_solve functions).
** gsl_sf_ellint_D  - remove useless argument n?
** Change default generator to Ranlxd (check seeding) or look at improved seeding for MT
** Vegas struct is too large and control variables should go at the beginning
** Remove use of long double internally, e.g. as an accumulator in
loops. It introduces variation between platforms which is undesirable.
It should be replaced with a preprocessor variable ACC_DOUBLE so
that the user can compile the library with the old long double
behavior if desired.
** Eliminate use of volatile where it has been used to force rounding
(integration/). It is better to write the code to avoid dependence on
rounding.
** Constant objects (like gsl_roots_fsolver_brent) ought to have
constant pointers (const gsl_roots_fsolver_type * const
gsl_roots_fsolver_brent)
** Make the return value EINVAL vs EDOM consistent for invalid
parameters. EDOM means a domain error (i.e. float or mathematically
undefined), EINVAL means invalid (i.e. zero length)
** Change return 0 to return GSL_SUCCESS, and return -1 to GSL_FAILURE
throughout, where appropriate. Similarly change any if(...) checks of
return values to use == GSL_SUCCESS, if they are checking for zero.
N.B. want to be careful about accidentally omitting error conditions
if using something like == GSL_FAILURE when function returns a
different error code.
** Make sure that all #defines are fully wrapped in ()'s, especially
the outermost layer which may have been missed. Everything should be
of the form #define foo(x) (....) so there is no possibility of bad
parsing.  Need a perl script to check this!
** Convert to BZR? (check GPG signing and integrity guarantees)

* Release "1.14" backwards compatible
** Add error code to all linalg svx type functions for singular matrix
** Merge integ-glfixed branch
** Fix bug  #25383 use GSL_ENOPROG instead of GSL_CONTINUE in lmiterate.c
** compare new erf function
** Export dwt.c/dwt_step
** Fix relative error in integration check

Re: GSL 2.0 roadmap

[Tuomo Keskitalo]

Hello,

thanks for the list. I was wondering that since there will be a new 
major version, this might be a good chance to modify some of the APIs / 
frameworks. Is this OK? For example, when I was implementing msadams and 
msbdf for ode-initval, I had to make some modifications to the 
algorithms in order to fit them into GSL ode-initval framework. It is 
clear that the stepper, control and evolve objects would sometimes 
benefit from mutual communication. Should ode-initval framework be 
changed so that this kind of communication is possible in future?

How long would it be before 2.0 is out, approximately? More or less than 
a year? It might be good idea to keep the 2.0 branch "unstable" for some 
time, before APIs are frozen (if changes to frameworks are OK).

On 08/21/2009 11:41 PM, Brian Gough wrote:

> * Changes for Release 2.0
> Break binary compatibility, but keep source compatibility.

Does "source compatibility" mean that user interfaces / APIs should not 
be changed?

> ** Convert to BZR? (check GPG signing and integrity guarantees)

Do you mean http://bazaar-vcs.org/Bzr ?
Would this mean to abandon git?

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: GSL 2.0 roadmap

[Robert G. Brown]

On Thu, 27 Aug 2009, Tuomo Keskitalo wrote:

> Hello,
>
> thanks for the list. I was wondering that since there will be a new major 
> version, this might be a good chance to modify some of the APIs / frameworks. 
> Is this OK? For example, when I was implementing msadams and msbdf for 
> ode-initval, I had to make some modifications to the algorithms in order to 
> fit them into GSL ode-initval framework. It is clear that the stepper, 
> control and evolve objects would sometimes benefit from mutual communication. 
> Should ode-initval framework be changed so that this kind of communication is 
> possible in future?
>
> How long would it be before 2.0 is out, approximately? More or less than a 
> year? It might be good idea to keep the 2.0 branch "unstable" for some time, 
> before APIs are frozen (if changes to frameworks are OK).

On the same note, a rather long time ago I developed a set of GSL-ish
"tensor" structures to enable the GSL to handle things beyond matrices
(which involved of course redefining slightly 1st and 2nd rank tensors,
or at least paralleling the matrix and the vector).  Tensors are,
obviously useful to all sorts of people doing computations in physics;
not only spatial tensors in e.g. relativistic theory computations but
also angular moment tensors that may not be efficiently rectangular.

Also, the GSL has lacked multidimensional quadrature; I ported one
algorithm from the C++ Hint package long ago to a GSLish form, but it
would be good at some point to consider adding an "extensible" set of
quadrature over d>1 spaces, where the "simple" solution of nested
integration is often horrendously inefficient and even inaccurate,
especially if the space involved is curved e.g. the surface of a sphere.

Finally, there are a number of possible additions to the roster of RNGs
out there, including relatively slow ones of cryptographic quality that
can serve purposes that even the best of the fast generators arguably
cannot, plus code that I wrote for dieharder that will automagically add
e.g.  /dev/random and /dev/urandom if they exist to the list.

I haven't worked with the first two for a while now, but I've still got
the code (of course) and would be happy to de-mothball it if there is
any interest in reconsidering the vector/matrix/tensor interface so it
is more extensible.  I think I implemented all the way out to 10th rank
tensors -- possibly overkill even for relativity, but I myself use at
least sixth rank in some of my code, and I'm aware of people who need
eighth rank.

    rgb

>
> On 08/21/2009 11:41 PM, Brian Gough wrote:
>
>> * Changes for Release 2.0
>> Break binary compatibility, but keep source compatibility.
>
> Does "source compatibility" mean that user interfaces / APIs should not be 
> changed?
>
>> ** Convert to BZR? (check GPG signing and integrity guarantees)
>
> Do you mean http://bazaar-vcs.org/Bzr ?
> Would this mean to abandon git?
>
> -- 
> Tuomo.Keskitalo@iki.fi
> http://iki.fi/tuomo.keskitalo
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap

[Jordi Burguet Castell]

Hello Robert,

On Thu, Aug 27, 2009 at 7:51 AM, Robert G. Brown<rgb@phy.duke.edu> wrote:
> On the same note, a rather long time ago I developed a set of GSL-ish
> "tensor" structures to enable the GSL to handle things beyond matrices
> (which involved of course redefining slightly 1st and 2nd rank tensors,
> or at least paralleling the matrix and the vector).  Tensors are,
> obviously useful to all sorts of people doing computations in physics;
> not only spatial tensors in e.g. relativistic theory computations but
> also angular moment tensors that may not be efficiently rectangular.
>
> [...]
>
> I haven't worked with the first two for a while now, but I've still got
> the code (of course) and would be happy to de-mothball it if there is
> any interest in reconsidering the vector/matrix/tensor interface so it
> is more extensible.  I think I implemented all the way out to 10th rank
> tensors -- possibly overkill even for relativity, but I myself use at
> least sixth rank in some of my code, and I'm aware of people who need
> eighth rank.

You may want to check the extensions "tensor" and "marray" that I also
developed a while ago:
http://www.gnu.org/software/gsl/#extensions

I don't know how much overlap they have with your ideas, the main
thing for them is to manipulate tensors of arbitrary rank and do the
common operations on them.

This may also interest you, quoting from the web: "If you want to add
a feature to GSL we recommend that you make it an extension first. We
will list it here so that people can try it out. Extensions can be
incorporated after they have been tested in real use (see "How to
help" for more information)."  So maybe you can consider making your
code an extension if it is not already?

Cheers,
Jordi

Re: GSL 2.0 roadmap

[Robert G. Brown]

On Thu, 27 Aug 2009, Jordi Burguet Castell wrote:

> You may want to check the extensions "tensor" and "marray" that I also
> developed a while ago:
> http://www.gnu.org/software/gsl/#extensions
>
> I don't know how much overlap they have with your ideas, the main
> thing for them is to manipulate tensors of arbitrary rank and do the
> common operations on them.

I'll try to look when I get a chance.  Unfortunately none of my current
projects are using tensors much these days, so it may take me some time
to get to it.  Basically my approach was to build a fairly
systematically creating ***...whatever pointers with arbitrary index
limits, creating a linear data object, and suitably packing the offsets
into the toplevel, castable tensor pointer.  The overall tensor object
was a struct with the requisite pieces and e.g. limit information
available inside the struct, with a create/destroy function that
allocated and packed and unpacked and freed respectively, but the
approach allowed one to directly access the tensor in array notation,
whatever[i][j][k]...  Not terribly subtle, in other words.  Your
extension is probably better thought out.

> This may also interest you, quoting from the web: "If you want to add
> a feature to GSL we recommend that you make it an extension first. We
> will list it here so that people can try it out. Extensions can be
> incorporated after they have been tested in real use (see "How to
> help" for more information)."  So maybe you can consider making your
> code an extension if it is not already?

Yes, I recall that, although it was a pretty long time ago when I wrote
these particular extensions and I'm not sure that this was yet standard
policy.  Maybe it was, I don't know; I do remember asking if anybody was
interested in them at that time and nobody suggested making them a
formal extension instead.  My dieharder program is listed in this sort
of way (including, I imagine, the extra GSL-compatible rngs and an
"add-rng" function that basically makes the GSL list extensible).  I
could probably turn at least that part into a formal extension pretty
easily, and it is arguably the most useful piece of what I've got.

Thanks!

    rgb

>
> Cheers,
> Jordi
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap

[Brian Gough]

At Thu, 27 Aug 2009 14:42:12 +0300,
Tuomo Keskitalo wrote:
> thanks for the list. I was wondering that since there will be a new
> major version, this might be a good chance to modify some of the
> APIs / frameworks. Is this OK?  For example, when I was implementing
> msadams and msbdf for ode-initval, I had to make some modifications
> to the algorithms in order to fit them into GSL ode-initval
> framework. It is clear that the stepper, control and evolve objects
> would sometimes benefit from mutual communication. Should
> ode-initval framework be changed so that this kind of communication
> is possible in future?

Yes, we could do that wherever there is a clear benefit.

> How long would it be before 2.0 is out, approximately? More or less than 
> a year? It might be good idea to keep the 2.0 branch "unstable" for some 
> time, before APIs are frozen (if changes to frameworks are OK).

Hopefully less than a year. I still try to make a release every six
months (more or less).

> Does "source compatibility" mean that user interfaces / APIs should not 
> be changed?

I think I meant to write "keep source compatibility where possible"
(i.e. for rearranging structs) .  For some changes, we will need to
break that.

> > ** Convert to BZR? (check GPG signing and integrity guarantees)
> 
> Do you mean http://bazaar-vcs.org/Bzr ?
> Would this mean to abandon git?

Bzr is now the standard GNU version control system, so it would mean
abandoning git, yes.  The timing of switching to git was unfortunate,
as the announcement about Bzr was not long after, but there wasn't any
way to know that before.  However, switching again isn't a high
priority.

-- 
Brian Gough

Re: GSL 2.0 roadmap

[Gerard Jungman]

Well, the discussion has started, and people are starting to step
out with their "2.0" ideas. The way I look at it, there are three
things that need to happen, in general terms:

  1) The backlog of 1.x-incompatible changes needs to be cleared.
     This corresponds to stuff on Brian's 2.0 todo list; probably
     lots of other stuff too. We need to make sure we overlook
     nothing.

  2) New modules can be considered for inclusion. 

  3) Changes to existing modules must be considered, where there are
     clear needs for improvements, in implementations, interfaces,
     or both.


Set (1) is fairly well-controlled, although it will take some
care to make sure everything is done right. I'm sure Brian is on
top of that. Set (2) is easy in a sense, since by definition
it involves well-encapsulated changes/additions.

Clearly, set (3) is the most difficult. Major issues here
are related to complex numbers, vector/matrix abstractions,
linear algebra, special functions, and FFTs. There are others,
but it will take some time to create a full assessment of the
existing modules. Detailed input from users will be needed.

However, we do need to decide on some plan for set (3), since
small changes in set (1) might affect or be affected by ideas
from (3). Some partial assessment of where we are is needed.

A GSL state-of-the-union is a large task. I have been thinking
about it, on and off, for months. In a followup message, I will
include my (incomplete) thoughts on this. These are almost all
critical thoughts, since I don't see any point in patting
ourselves on the back.

Feel free to blast away at it.

--
Gerard Jungman

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

[-- Attachment #1: Type: text/plain, Size: 1 bytes --]



[-- Attachment #2: gsl-2.0-outline-2009.08.27.txt --]
[-- Type: text/plain, Size: 15443 bytes --]

Requirements (a very incomplete list)
-------------------------------------

* A controlled shift in the C interfaces. Too much interface change
  may scare off 1.x users. But there must be change...

* Better integration with existing toolsets.
  - vector / matrix models
  - linear algebra (!)

* Better overall organization, expression of module dependencies, etc.

* GSL should be a technology leader. That is part of the point of the
  GPL licensing. If GSL were just another of many essentially equivalent
  libraries, then it would make more sense for it to be LGPL or even BSD.
  To justify GPL licensing, which coerces people into our open-source
  world, we better have something unique and excellent.



Taking Stock of 1.x
===================


** Overall Design

 * Split into "independent sublibraries".

   The stated design goal of independent sublibraries was abandoned
   in all but name quite early in the project. Several obvious problems
   with a naive interpretation of "independent" occur:

    - Some sublibraries depend on other "foundational" sublibraries;
      such dependencies are natural and desirable. But with no effective
      explicit way to express these dependencies, the overall organization
      degenerates to a monolith.

    - The packaging environment does not _easily_ support independence.
      This is related to the lack of an effective expression of
      dependencies among sublibraries, as noted above. It is also
      related to practical problems in supporting a hierarchical
      build with the autotools setup. Witness the apparent need to
      create the gsl/ directory of links to header files as the
      first step in the build hackery.

   Consequences: Various practical problems, such as the organization
                 of header files, the grotesque link lines (a real problem
		 at one point), and general difficulty in comprehension of
		 the structure of the project.

   Requirement: A natural, effective, and explicit way to express
                dependencies, which can be exported to the build
		and packaging systems

 * One Language Only

   The benefits of a one-language only design are described in the
   design document. They are fairly obvious and depend on the fact
   that C is the "universal" system language.

   Nevertheless, there are large costs associated to this choice.
   Amongst these are the following.

    - Important "legacy" tools which are implemented in other languages
      are unavailable to GSL. This creates a major deficiency in areas
      such as linear algebra. Some holes were partially plugged with
      native GSL implementations, but this can never be an acceptable
      solution, either from a software design standpoint or from an
      end-user standpoint. From a software-design standpoint, GSL fails
      to gain from the maturity and the ongoing development in these
      external tools. From an end-user standpoint, the lack of performance
      has caused many users to abandon these aspects of GSL. Further
      discussion of these points occurs in specialized sections below,
      as appropriate.

    - Many new developments in methodology are occuring outside C, notably
      in C++. As long as such developments were essentially experimental and
      unfocused, they could be ignored; one could argue that this was indeed
      the case a decade ago, when GSL was initiated. But this can no longer
      be argued. The performance of these new tools exceeds (often by a
      large factor) the performance of the GSL tools; some are
      well-established, readily available, and of clear utility.
      Again, many end-users have abandoned the corresponding areas of GSL.

    - When a user abandons some core area of GSL for other implementations,
      they tend to abandon other aspects as well, since GSL interfaces are not
      always friendly to other data models. In many cases it would be
      impossible to make it so, due to insuperable inter-language barriers.
      GSL should not continue to ignore these developments. At the very least
      it should allow for the existence of these external developments, by
      paying very close attention to inter-language issues.

    These issues are related to the design goal of "naturalness", stated in
    the design document. Quote:
      "If there is something which is unnatural in C and has to be simulated
       then we avoid using it."
    It may be coming to pass that some/many useful methods and designs for a
    numerical library are now "unnatural in C" by onstruction, since the performant
    and well-designed tools are no longer being created in C. For example, it
    may be the case that problem domains such linear algebra are no longer
    "natural in C". Whither GSL in such a world?


 * Error Handling

   GSL implements a very rudimentary type of error-handling, which is
   quite appropriate for small projects. The "register handler, fail-by-default"
   model works for a small library with a "low stack depth" usage model. By this,
   I mean that library functions are called by clients, but not by other library
   functions; in such a case, it is easy for the client to interpret the errors.
   But GSL does not conform to this "low stack depth" usage model, and the
   fail-by-default model can be confusing for clients, when the failure
   occurs at depths that they cannot control.

   Although it may not seem like a pressing issue, this monolithic
   error-handling design seems to be another symptom of the lack of
   explicit models for interdependence in GSL. It should be revisited.



** Build System

   The autotools build system can be an impediment to overall design.
   As discussed above, the GSL build architecture, based on a typical
   autotools single-project setup, could not be coerced into supporting
   notions of dependence/independence for sublibraries. Rather, these
   dependencies exist in spite of the build system.



** Complex Numbers

   GSL implements complex numbers. This functionality may be superfluous at
   this time, depending on the status of C99 features in the compilers of
   interest.

   Beyond this observation, there are some thorny issues with complex numbers,
   when one considers inter-language issues. As argued in several places here,
   GSL cannot afford to ignore these issues. Users are often turning to basic
   functionality implemented in other languages (linear algebra!), and
   inter-operability at the level of data interchange is an absolute
   necessity. An array of complex values computed in GSL should be
   appropriate for passage through an inter-language interface.

   The two languages of most concern are fortran and c++. In fortran, complex
   numbers are a language feature, and in c++ they are a library feature.
   In both cases, the representation is, in principle, implementation-
   dependent. These problems have been addressed by other projects;
   effective solutions exist for the interface to fortran. It seems
   likely that the c-c++ interface problem is solvable, especially
   within the context of a consistent build with a single compiler.
   Mixing compilers may be more difficult, but should be solvable
   in principle.



** Vector and Matrix Data Structures [ vector/, matrix/, block/ ]

  The basic features of the GSL vector and matrix data structures
  are reasonable, as an implementation of dense storage for basic
  data types.

  The important notion of slicing is (partially) implemented in GSL
  in terms of the "view" concept. One can construct submatrices as
  views of given matrices, change the stride of vector data by
  creating vector views, etc. But there are clear flaws in the
  design. The design does not express the obvious idea that
  a "view" is itself a "thing", simply because the view classes
  do not have an inheritance relationship to the main classes.

  In fact, there is a problem with the functoriality of the design,
  because there is an obvious logical sense in which a matrix is
  itself a view type, using a view notion corresponding to the
  default strides, but the design seems to "point the other way".

  There are also significant usability issues with the simple aspects
  of the interfaces. The get() and set() functions have been deplored
  by many users. This syntax leads to unwieldy user code to accomplish
  the simplest tasks. The unwieldy nature of these interfaces often
  causes the user to introduce otherwise unnecessary temporaries.
  Although the compiler can remove such temporaries, the user cannot,
  and they end up occupying valuable _intellectual_ real estate.

  These problems stem fundamentally from the fact that the view
  classes were added after the main classes, and the design of the
  main classes was not modified as logic would have dictated. This
  is a failure of design.

  Executive Summary: I was there when this stuff was born,
  and I don't even understand it. I have to stare at the header
  files whenever I try to do anything with vectors.



** Basic Vector Matrix Operations [ blas/ ]

   The existence of a C interface standard for BLAS functionality
   greatly simplified the GSL apporach to BLAS. In GSL, all BLAS
   function usage adheres to the cblas interface standard, meaning
   that any cblas conformant binary interface is acceptable; any
   library which exports a cblas binary interface can be linked
   to a GSL application. For example, the ATLAS library provides
   a full cblas conformant implementation and can be linked to
   any GSL application.

   GSL also provides a native implementation of the cblas interfaces, based
   on the reference implmentation for the standard, which was a draft
   standard at the time of GSL 1.0 release. This is where GSL blas
   functionality begins to fall down. The GSL cblas implementation was felt
   to be a necessary evil, simply because not all installations could be
   expected to have a cblas conformant library installed. Unfortunately,
   many users are unaware of the underlying architecture and end up using
   the native GSL implementation by default, resulting in poor performance,
   which reflects badly on GSL as a whole.

   This is a natural place where an inter-language approach would have been
   useful. The standard fortran BLAS implementation is readily available and
   of adequate performance for almost all users. For example, in the current
   era, it is readily available as an rpm package. It would have been natural
   at that time to use cblas wrappers over the underlying fortran BLAS
   (wrappers perhaps also available as part of the draft cblas standard). For
   completeness, this may have required shipping the fortran source as a
   necessary part of GSL; this idea was rejected because it violated the "One
   Language" design rule. As a consequence, the developers responsible for the
   actual implementation (mainly myself!) were forced to transcribe the cblas
   reference implementation into GSL, resulting in about 8000 lines of needless
   code, violating every reasonable notion of software reuse.

   In a case like cblas, where an interface standard exists, the link-time choice
   method currently implemented in GSL seems like the correct solution. 

   When faced with the choice of using a decent external implementation in a 
   language other than C and re-implementing (badly) in C, "One Language"
   religious issues should not take precedence in the decision process.


 * sublibs that use blas
   - eigen/
   - linalg/
   - multifit/
   - multimin/



** Linear Algebra [ linalg/ ]

   The defacto standard for linear algebra functionality is LAPACK. LAPACK
   exists as (and is essentially defined by) a fortran implementation which
   is commonly available, on essentially all platforms. Like blas, I can
   just 'yum install' it.

   Partly because of the lack of a C standard interface, and partly because
   of "One Language" religious beliefs, this defacto standard implementation
   was not available for the design of GSL. Rather, GSL relies on a native
   implementation of some (smallish) subset of lapack functionality.

   The resulting lack of performance has driven many users away from
   GSL. The often heard advice on the GSL discussion list is that
   "serious users" should "use LAPACK instead". But which users
   of GSL consider themselves to be "non-serious"? The notion is
   absurd. It is not users which lack seriousness, but the
   GSL implementation.

   Linear algebra is a foundational tool in numerical computing.
   Users need it; other modules in GSL need it. We need a
   serious solution for this mess.

 * sublibs that use linear algebra
   - eigen/
   - interpolation/
   - multifit/
   - multiroots/
   - ode-initval/



** Random Number Generation

   The GSL random number generators represent one of the few foundational
   aspects of the library which can be considered successful. RNGs naturally
   lend themselves to a kind of shallow object-oriented design which is
   very natural in C. The design allows for easy extension, and a user
   who understands one of the rng types essentially understands them all.

   I believe this success follows from the relative simplicity of the
   problem domain, as far as it bears on the interface design.



** Special Functions

   The special function sublibrary in GSL suffers mainly from a lack
   of consistency and coherence in the implementations. This stems
   from the decision made early in the project that coverage was
   very important; if users had to leave GSL (and use that _other_ library)
   to evaluate functions, they would likely leave GSL aside entirely.
   Correctness was not sacrificed in principle, but in practice the
   implementations span the spectrum, from ironclad to somewhat-fishy.

   As a kind of apology for this state of affairs, the functions attempt
   to estimate errors, using heuristic and sometimes ill-defined methods.
   This turns out to be a poor apology, since it tends to gum-up the
   works for the whole sublibrary, eating performance and occupying
   a large piece of intellectual real estate.

   The error estimation code must, at the least, be factored out. More
   to the point, is should likely be discarded in the main line of
   development. Other notions of error control should be investigated.

   Similar to the notion of sublibrary dependence discussed elsewhere,
   the special functions should be hierarchically organized. The
   dependencies should be made clear, and foundational levels of
   ironclad functions should be made explicit. Such ironclad
   functions should occupy the same place in the users mind as
   platform-native implementations of sin(x); they should be
   beyond question for daily use.

   Other functions, which suffer from implementation problems
   or are simply too complicated to guarantee the same level
   of correctness should be explicitly identified, as part
   of the design (not just the documentation!).


** FFT

   The GSL native fft implementations suffer in the same way that
   the GSL native linear algebra implementations suffer. Better
   solutions are available, from almost all points of view, and
   the GSL implementations therefore detract from the GSL
   cachet as a whole. "Serious users should use FFTW" is not
   a good tagline for the sublibrary; see the discussion of
   linear algebra above.

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

Interesting...  it may take me a while to read this properly and reply
as I am on the road now.  I guess the main point I would add is that
we have to operate under constraints, compared to LAPACK etc we don't
have people funded to work on these things the way they do.  So there
is a question of whether we want to try to change that aspect, or just
make technical changes. 

-- 
Brian 

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Fri, 2009-08-28 at 14:48 +0100, Brian Gough wrote:

> I guess the main point I would add is that
> we have to operate under constraints, compared to LAPACK etc we don't
> have people funded to work on these things the way they do.  So there
> is a question of whether we want to try to change that aspect, or just
> make technical changes. 

This is an important point. But you seem to have it upside down.
It's precisely because we don't have these resources that we need
to use what other people have done, i.e. to _use_ things like
lapack, and make it easy for users to use lapack and gsl
simultaneously, without head scratching.

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Thu, 27 Aug 2009 17:15:39 -0600,
Gerard Jungman wrote:
> * Better overall organization, expression of module dependencies, etc.

I would be happy to remove the goal of independent sublibraries from
the design document--it doesn't reflect reality. 

I am not clear what benefit it would provide to users now anyway--is
there actually a problem with "one big library"?  The original
motivation was to minimise memory usage and allow code to be
extracted.  

Memory for shared libraries is not really an issue these days--we can
also leave it up to the linker to remove unnecessary functions if that
is a problem.  The idea of people extracting parts of the library
doesn't seem important either -- it is much easier to install the
whole thing.

If you want to avoid the symbolic linking of header files we could
just move them into the gsl/ directory instead of keeping them in the
individual subdirectories. (While the symlinks are ugly I don't think
they have ever actually caused a problem for anyone.)

-- 
Brian

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Thu, 2009-09-03 at 19:19 +0200, Brian Gough wrote:
> At Thu, 27 Aug 2009 17:15:39 -0600,
> Gerard Jungman wrote:
> > * Better overall organization, expression of module dependencies, etc.
> 
> I would be happy to remove the goal of independent sublibraries from
> the design document--it doesn't reflect reality. 

There are really several (partially) independent issues, all
touched by the "Better organization" point. They relate to
the build system, module dependencies, and the overall user
experience.


> I am not clear what benefit it would provide to users now anyway--is
> there actually a problem with "one big library"?

Yes, there is. As soon as GSL learns to talk with the rest of the
world, dependency issues will arise. Sub-libraries should be
insulated. Think of how something like a dependency on lapack
might propagate through the system.

There are related problems, related to dependencies and
the nature of the build, which can be addressed by
compartmentalizing the installation. This is helpful
to package maintainers, whether they are preparing packages
for a major distribution or they are local ad-hoc
packagers, down to the level of one user.

Simply put, the build and installation should be more modular.


> If you want to avoid the symbolic linking of header files we could
> just move them into the gsl/ directory instead of keeping them in the
> individual subdirectories. (While the symlinks are ugly I don't think
> they have ever actually caused a problem for anyone.)

This is the upside down solution. The right way is simple; just change
the includes to be hierarchical, so that client code would look like

  #include <gsl/interpolation/gsl_spline.h>

We talked about this 12 years ago, but we couldn't figure out how
to get autotools to do this properly, for build and/or installation.
Maybe somebody knows how. Maybe we should discard autotools.
That should certainly be discussed.

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Thu, 27 Aug 2009 17:15:39 -0600,
Gerard Jungman wrote:
>   There are also significant usability issues with the simple aspects
>   of the interfaces. The get() and set() functions have been deplored
>   by many users. 

Is there an alternative in C?  I can only think of providing macros
for expressions like *gsl_vector_ptr(v,i) = x if one prefers to assign
to an l-value.

-- 
Brian Gough

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Mon, 2009-09-07 at 16:07 +0100, Brian Gough wrote:
> At Thu, 27 Aug 2009 17:15:39 -0600,
> Gerard Jungman wrote:
> >   There are also significant usability issues with the simple aspects
> >   of the interfaces. The get() and set() functions have been deplored
> >   by many users. 
> 
> Is there an alternative in C?  I can only think of providing macros
> for expressions like *gsl_vector_ptr(v,i) = x if one prefers to assign
> to an l-value.

Yes and no.

For example, in 9/10 cases the underlying data is just a contiguous
block of double (with unit stride). In that case I can always just
extract the pointer and apply [].

Ok, you say, this is very dangerous because somebody could pull the
rug out from under you, or you might forget about striding, etc.
True. But the GSL interface should not make it harder to do things
in the simplest 9/10 cases. There should be a way for the data
pointer (or block) to coexist with the meta-info for layout,
so that both can be used as needed.

It smells like a view, and indeed it is. But it is a notion of
view which is more central to the design than the current one.
Everything should be a view.

Again, more details will follow, if I get around to typing
what is in my head. But we're all experienced adults here;
isn't it clear that there must be a better way, independent
of what I end up typing?

--
G. Jungman

Re: new double precision data structure?

[Tuomo Keskitalo]

[-- Attachment #1: Type: text/plain, Size: 674 bytes --]

On 09/16/2009 03:49 AM, Gerard Jungman wrote:

> It smells like a view, and indeed it is. But it is a notion of
> view which is more central to the design than the current one.
> Everything should be a view.

I think we really need to sort the data structure out now.

I'm not sure what you have in mind, but to get some discussion going on, 
I've attached an (incomplete and untested) idea for a general 
n-dimensional double precision data structure. Is this anything like 
what you pictured?

I'd like to hear people's comments on this. Would something like this be 
of use? Feel free to give your full critique.

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

[-- Attachment #2: newblock.c --]
[-- Type: text/x-csrc, Size: 4355 bytes --]

/* newblock.c - An idea for n-dimensional double precision data
   structure. This is pseudocode to describe the idea, does not
   compile.
*/

/* Structure for n-dimensional double precision data storage */

typedef struct 
{
  size_t dim;     /* data structure dimension */
  size_t * sizes; /* dimensions of system, e.g. {3,4,5} for real 3*4*5 tensor */
  double * data;  /* continuous data storage */
  size_t layout;  /* data layout type - LAST_INDEX_CONTINUOUS =
		     "column major" for matrices, or
		     FIRST_INDEX_CONTINUOUS = "row major" (default)
		     for matrices. Note: It might be possible to
		     extend this to support other storage schemes,
		     e.g. diagonal/banded storage
		     BANDED_INDEX_CONTINUOUS. */
} base_block;

/* Structure for a generic n-dimensional view object to access base
   block data. Access to data storage is done via this view
   object. There can be multiple views referring to the same data
   storage.
*/

typedef struct 
{
  size_t * firstcell; /* indices to first access cell (minimum indices) */
  size_t * lastcell; /* indices to last access cell (maximum indices) */
  size_t * stride; /* strides for each dimension for the view */ 
  base_block * block;
  double * data; /* pointer to block->data */
  size_t owner; /* indicator for view owning the data */
} base_view;

/* vector object based on base block and view */

typedef struct 
{
  base_view * view;
} vector;

/* matrix object based on base block and view */

typedef struct 
{
  base_view * view;
} matrix;

matrix * matrix_alloc (size_t n1, size_t n2) 
{
  /* Allocates matrix object (both view and storage objects) in
     row-major layout. */
  // 1. allocate base_block with dim=2, sizes={n1,n2}
  // 2. allocate base_view with firstcell={1,1}, 
  //    lastcell={n1,n2}, stride={0,0} (=whole matrix)
  // 3. set view owner=TRUE
  // 4. return base_view;
}

matrix * matrix_alloc_with_layout (size_t n1, size_t n2, size_t layout) 
{
  /* Allocates matrix object (both view and storage objects) in the
     specified layout. */
}

int vector_alloc_view_of_matrix (vector * v, matrix * m, 
				 size_t * firstcell, size_t * lastcell, 
				 size_t * stride) 
{
  /* Creates a new vector view of the matrix m */
  // 1. allocate base_view with given firstcell, lastcell and strides
  // 2. associate view with the storage
  // 3. set owner=FALSE
  // 4. return base_view;
}

/* Example code */

int main (void) 
{
  matrix * a;
  matrix * b;
  vector * v;
  int retval;

  /* allocate matrix with 3 rows and 4 columns */
  a = matrix_alloc (3, 4);

  /* User can change the layout between row- and column-major at any
     point. After change to column-major, those GSL functions that
     require row-major storage would raise error if the matrix is
     passed to them.
  */
  retval = matrix_set_layout (a, LAST_INDEX_CONTINUOUS); 

  {
    double temp[] = { 0, 10, 20, 01, 11, 21, 02, 12, 22, 03, 13, 23 };
    size_t templen = 12;
    retval = matrix_set_doublevec (a, temp, templen);
    /* This would initialize matrix
       a = [00 01 02 03
            10 11 12 13
            20 21 22 23] */
  }

  /* Set up a vector view of the diagonal */
  const size_t b0[] = { 1, 1 };
  const size_t e0[] = { 3, 3 };
  const size_t s0[] = { 1, 1 };
  retval = vector_alloc_view_of_matrix (v, a, &b0, &e0, &s0);
  
  /* Access cells in matrix through the vector */
  double val;
  val = vector_get(v, 1); // returns 00 from matrix a
  val = vector_get(v, 2); // returns 11 from matrix a
  val = vector_get(v, 3); // returns 22 from matrix a
  retval = vector_set(v, 2, 55); // sets value 55 in place of 11 to matrix a

  size_t l;
  l = vector_length(v); // returns 3;

  double diag[3];
  retval = vector_get_doublevec (&diag, v); // copies diagonal element values to diag
}

/* complex number matrix (one possible implementation) based on base
   block and view */

typedef struct 
{
  base_view * view;
} complex_matrix;

complex_matrix * complex_matrix_alloc (size_t n1, size_t n2) {
  /* Allocates matrix object (both view and storage objects) */
  // 1. allocate base_block with dim=3, sizes={n1,n2,2}, so that 
  //    the last dimension spans the real and imaginary parts
  // 2. allocate base_view with firstcell={1,1,1} and 
  //    lastcell={n1,n2,2} (=whole matrix)
  // 3. set view owner=TRUE
  // 4. return base_view;
}

Re: new double precision data structure?

[James Bergstra]

I really like the idea of a GSL n-dimensional array, because it would
make it much easier to interoperate with numpy, scipy, and the other
tools which are developed in that community.

Compared to numpy's ndarray structure, there are only a few
differences with your proposal.

Firstly, the ndarray untyped.  The data is in a void * pointer or a
char * pointer or something, and there is an extra enum-valued field
that indicates what sort of elements make up the data.  For example, 0
might mean int8, 1 might mean uint8, 2 might mean int32, 10 might mean
float32, 11 float64, 12 complex64, and so on.  There is support for
non-native elements like structs and objects too.  We wouldn't need
all of that flexibility in the GSL, but the technique might be useful
for allocating and passing around tensors of integers, floats, complex
numbers.

Secondly, there is no notion of a 'layout' in the ndarray, so I wonder
if it is necessary here?  When you explicitly store the strides for
traversing the tensor in each dimension... what is left for the
'layout' to specify?

Thirdly, the dimensions are logical things, not physical ones.  So
perhaps they belong in the view rather than the base block?  I think
the simple { size_t n_allocated; char * buf } structure is sufficient
for the base block.

Fourthly, if you take the dimensions out of the base block, then they
need to into the view as well... so the view would wind up looking
something like:

{
int ndim;
size_t * dims; /* the size of the tensor in each of `ndim` dimensions */
int * strides;  /*negative strides are very useful */
gsl_type_t type; /* something like gsl_float32, gsl_float64,
gsl_complex128, etc. */
void * data;  /* pointer into the underlying block, not necessarily
the beginning */
base_block * base;
int owner; /* True -> free the base with ourselves */
}

The biggest advantage of moving the dimensions from the base_block to
the view that it allows in-place reshaping and transposing of views,
which would be awkward otherwise.

Macros / functions can make it convenient to get a properly-casted
pointer to the underlying data (returning NULL when the type being
asked for doesn't match the type of the data).  For example
gsl_view_data_float32(view).

James Bergstra
-- 
http://www-etud.iro.umontreal.ca/~bergstrj


On Sun, Sep 27, 2009 at 4:03 AM, Tuomo Keskitalo <Tuomo.Keskitalo@iki.fi> wrote:
> On 09/16/2009 03:49 AM, Gerard Jungman wrote:
>
>> It smells like a view, and indeed it is. But it is a notion of
>> view which is more central to the design than the current one.
>> Everything should be a view.
>
> I think we really need to sort the data structure out now.
>
> I'm not sure what you have in mind, but to get some discussion going on,
> I've attached an (incomplete and untested) idea for a general n-dimensional
> double precision data structure. Is this anything like what you pictured?
>
> I'd like to hear people's comments on this. Would something like this be of
> use? Feel free to give your full critique.
>
> --
> Tuomo.Keskitalo@iki.fi
> http://iki.fi/tuomo.keskitalo
>

Re: new double precision data structure?

[Tuomo Keskitalo]

Hello,

thanks, good comments!

On 09/27/2009 07:34 PM, James Bergstra wrote:

> Firstly, the ndarray untyped.  The data is in a void * pointer or a
> char * pointer or something, and there is an extra enum-valued field
> that indicates what sort of elements make up the data.  For example, 0
> might mean int8, 1 might mean uint8, 2 might mean int32, 10 might mean
> float32, 11 float64, 12 complex64, and so on.  There is support for

This implies that GSL should support other elementary types than double. 
What do people think about this?

I rather like the notion. Maybe we can't provide all GSL functionality 
to all those data types, but it would certainly be nice to have this 
option on the background, even if double precision would be the only one 
fully supported.

> Secondly, there is no notion of a 'layout' in the ndarray, so I wonder
> if it is necessary here?  When you explicitly store the strides for
> traversing the tensor in each dimension... what is left for the
> 'layout' to specify?

The fact whether e.g. matrix data is stored in row or column major 
format. For example, if m is 2*2 matrix stored in d, and if d[0] is the 
(1,1) cell in matrix, then d[1] might be (2,1) (=column rajor storage, I 
think) or (1,2) (=row major storage), depending on your programming 
language/library.

> Thirdly, the dimensions are logical things, not physical ones.  So
> perhaps they belong in the view rather than the base block?  I think
> the simple { size_t n_allocated; char * buf } structure is sufficient
> for the base block.

I was thinking that the actual dimensions would be stored in the block, 
and the view would define an allowed subspace. For example, a view for a 
matrix could be a matrix or a vector. So, the geometry of block level 
would be frozen, and the views would provide alternative ways to access 
the data in it in different geometries. Is this too restrictive?

> The biggest advantage of moving the dimensions from the base_block to
> the view that it allows in-place reshaping and transposing of views,
> which would be awkward otherwise.

I did not get this, could you please give an example? If there is a 
problem, then dimensions can go to the view as you suggest.

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: new double precision data structure?

[Brian Gough]

At Mon, 28 Sep 2009 18:47:50 +0300,
Tuomo Keskitalo wrote:
> On 09/27/2009 07:34 PM, James Bergstra wrote:
> > Firstly, the ndarray untyped.  The data is in a void * pointer or a
> > char * pointer or something, and there is an extra enum-valued field
> > that indicates what sort of elements make up the data.  For example, 0
> > might mean int8, 1 might mean uint8, 2 might mean int32, 10 might mean
> > float32, 11 float64, 12 complex64, and so on.  There is support for
> 
> This implies that GSL should support other elementary types than double. 
> What do people think about this?

I think the types in multidimensional arrays ought to follow the
existing arrangement and naming convention as for vector and matrices.

Re: new double precision data structure?

[Gerard Jungman]

On Sun, 2009-09-27 at 11:03 +0300, Tuomo Keskitalo wrote:
> 
> I'm not sure what you have in mind, but to get some discussion going on, 
> I've attached an (incomplete and untested) idea for a general 
> n-dimensional double precision data structure. Is this anything like 
> what you pictured?

I started working on something over the weekend; hope to have
a working prototype (at least for double) in a week or so. I
will have a look at your thing now and see what i can steal
from it...

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Thu, 27 Aug 2009 17:15:39 -0600,
Gerard Jungman wrote:
>   The important notion of slicing is (partially) implemented in GSL
>   in terms of the "view" concept. One can construct submatrices as
>   views of given matrices, change the stride of vector data by
>   creating vector views, etc. But there are clear flaws in the
>   design. The design does not express the obvious idea that
>   a "view" is itself a "thing", simply because the view classes
>   do not have an inheritance relationship to the main classes.

I agree that the scheme is not as elegant as it could be in other
languages.  The view types are forced by the nature of const in C --
it's not possible to place the views and vectors/matrices on an equal
footing and preserve constness, unfortunately.  If there is a way
round that, I'm not aware of it.

Originally a view was essentially a vector with another name, but
there was no barrier to writing expressions which discarded constness
without any complaint from the compiler.  To prevent that in C, the
type had to be "wrapped" in a struct which is why one has to write
&view.vector or &view.matrix to access it.

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Mon, 2009-09-07 at 16:06 +0100, Brian Gough wrote:
> At Thu, 27 Aug 2009 17:15:39 -0600,
> Gerard Jungman wrote:
> >   The important notion of slicing is (partially) implemented in GSL
> >   in terms of the "view" concept. One can construct submatrices as
> >   views of given matrices, change the stride of vector data by
> >   creating vector views, etc. But there are clear flaws in the
> >   design. The design does not express the obvious idea that
> >   a "view" is itself a "thing", simply because the view classes
> >   do not have an inheritance relationship to the main classes.
> 
> I agree that the scheme is not as elegant as it could be in other
> languages.  The view types are forced by the nature of const in C --
> it's not possible to place the views and vectors/matrices on an equal
> footing and preserve constness, unfortunately.  If there is a way
> round that, I'm not aware of it.

> Originally a view was essentially a vector with another name, but
> there was no barrier to writing expressions which discarded constness
> without any complaint from the compiler.  To prevent that in C, the
> type had to be "wrapped" in a struct which is why one has to write
> &view.vector or &view.matrix to access it.


I really don't understand this. In my head, I can see a solution
which has nothing to do with const-ness issues. I think it would
just work. I could type it in and we could look at it. It definitely
involves changing the way vector and matrix are done, but I don't
think it would be a big hairy deal.


On Mon, 2009-09-07 at 14:21 -0400, Robert G. Brown wrote:
>
> The biggest issue would/will probably be rationalizing the views of
> vector and matrix so they are sufficiently portable and easy to e.g.
> pass in and out of ODE solvers and everything else consistently.

This is right to the point. It's exactly what I'm getting at.
There is just something brain-damaged about the way it is done now.

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Robert G. Brown]

On Tue, 15 Sep 2009, Gerard Jungman wrote:

> On Mon, 2009-09-07 at 16:06 +0100, Brian Gough wrote:
>> At Thu, 27 Aug 2009 17:15:39 -0600,
>> Gerard Jungman wrote:
>>>   The important notion of slicing is (partially) implemented in GSL
>>>   in terms of the "view" concept. One can construct submatrices as
>>>   views of given matrices, change the stride of vector data by
>>>   creating vector views, etc. But there are clear flaws in the
>>>   design. The design does not express the obvious idea that
>>>   a "view" is itself a "thing", simply because the view classes
>>>   do not have an inheritance relationship to the main classes.
>>
>> I agree that the scheme is not as elegant as it could be in other
>> languages.  The view types are forced by the nature of const in C --
>> it's not possible to place the views and vectors/matrices on an equal
>> footing and preserve constness, unfortunately.  If there is a way
>> round that, I'm not aware of it.
>
>> Originally a view was essentially a vector with another name, but
>> there was no barrier to writing expressions which discarded constness
>> without any complaint from the compiler.  To prevent that in C, the
>> type had to be "wrapped" in a struct which is why one has to write
>> &view.vector or &view.matrix to access it.
>
>
> I really don't understand this. In my head, I can see a solution
> which has nothing to do with const-ness issues. I think it would
> just work. I could type it in and we could look at it. It definitely
> involves changing the way vector and matrix are done, but I don't
> think it would be a big hairy deal.
>
>
> On Mon, 2009-09-07 at 14:21 -0400, Robert G. Brown wrote:
>>
>> The biggest issue would/will probably be rationalizing the views of
>> vector and matrix so they are sufficiently portable and easy to e.g.
>> pass in and out of ODE solvers and everything else consistently.
>
> This is right to the point. It's exactly what I'm getting at.
> There is just something brain-damaged about the way it is done now.

Bless you.  And remember, it would be really nice if there was a certain
scalability in dimension.  The world doesn't stop at *a, **b -- there is
***c, ****d, etc where a number of sciences need at least through six
dimensions and back when we were kicking this around on list three or
four years ago somebody (doing relativity?) said that they use eight or
nine.  Tensors, especially tensors that don't "have" to be square, are
very useful.

In my own code, I have used at least ****d tensors that were allocated
as a single data block and the indices repacked into the tensor and then
cast to a type.  You can bass the data block to an ODE solver as a
vector, but you can write the actual deriv routine using the tensor
pointer, which means that it is actually intelligible and debuggable
without four layers of index arithmetic wrapped up on some sort of macro
or horrorshow code.

I don't remember if I mentioned it before, but multidimensional
integration (on spaces that may or may not be flat of dimension greater
than or equal to two) would be a lovely addition to the next release as
well.  There's C++ code out there (Hint) and I did a port of one of the
Hint routines at one time, but e.g. integration over a sphere of
functions with Y_lm's in them is another extremely common problem in
math and physics, and spherical domains in spherical coordinates do not
map well into nested one-D quadrature routines...

    rgb

>
> --
> G. Jungman
>
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Tue, 15 Sep 2009 18:49:07 -0600,
Gerard Jungman wrote:
> I really don't understand this. In my head, I can see a solution
> which has nothing to do with const-ness issues. I think it would
> just work. I could type it in and we could look at it. It definitely
> involves changing the way vector and matrix are done, but I don't
> think it would be a big hairy deal.

If you can make a simple concrete implementation of how it could work
that would help a lot.

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Thu, 27 Aug 2009 17:15:39 -0600,
Gerard Jungman wrote:
>    GSL implements a very rudimentary type of error-handling, which
>    is quite appropriate for small projects. The "register handler,
>    fail-by-default" model works for a small library with a "low
>    stack depth" usage model. By this, I mean that library functions
>    are called by clients, but not by other library functions; in
>    such a case, it is easy for the client to interpret the errors.
>    But GSL does not conform to this "low stack depth" usage model,
>    and the fail-by-default model can be confusing for clients, when
>    the failure occurs at depths that they cannot control.

The issue here is threading -- otherwise we could use a global
variable to control the behaviour.  

Either we have to 

- pass an extra argument to all(?) functions for error control (seems
  a large impact for a relatively small problem)

- provide separate _impl non-aborting versions for functions called
  internally, similar to the special functions (simplest option).

- use thread local storage for a global error control variable (not
  100% portable).

-- 
Brian Gough

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

regarding error handling:

On Mon, 2009-09-07 at 14:57 +0100, Brian Gough wrote:
> 
> The issue here is threading -- otherwise we could use a global
> variable to control the behaviour.  

> Either we have to 
> 
> - pass an extra argument to all(?) functions for error control (seems
>   a large impact for a relatively small problem)
> 
> - provide separate _impl non-aborting versions for functions called
>   internally, similar to the special functions (simplest option).

A combination of these seems right. We also have the usual problem
of error codes for functions with natural prototypes. We should
finally sort that out and do something coherent across the library.
The _e and natural versions of the specfunc functions are just one
example of that. Any ideas are welcome. We should think about this
now, since it affects the interfaces across the library.


> - use thread local storage for a global error control variable (not
>   100% portable).

I would like to know how non-portable this is. I mean, I understand it
is different for every platform, but if everybody has one then it
only becomes an issue of configuring the variability.

Roughly speaking, I would think that every thread implementation must
have a notion of thread-local storage, so we're never required to
provide something that doesn't exist...

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Tue, 15 Sep 2009 18:48:32 -0600,
Gerard Jungman wrote:
> I would like to know how non-portable this is. I mean, I understand it
> is different for every platform, but if everybody has one then it
> only becomes an issue of configuring the variability.
> 
> Roughly speaking, I would think that every thread implementation must
> have a notion of thread-local storage, so we're never required to
> provide something that doesn't exist...

I don't know enough about it to really say, but I haven't seen any
libraries following that approach -- returning an error code is the
usual way.

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

>    As a kind of apology for this state of affairs, the functions attempt
>    to estimate errors, using heuristic and sometimes ill-defined methods.
>    This turns out to be a poor apology, since it tends to gum-up the
>    works for the whole sublibrary, eating performance and occupying
>    a large piece of intellectual real estate.
> 
>    The error estimation code must, at the least, be factored out. More
>    to the point, is should likely be discarded in the main line of
>    development. Other notions of error control should be investigated.

I have found the error estimates to be pretty useful when debugging,
and they are only a small extra cost.

Some of the hypergeometric functions are difficult to understand in
places.  Here it would be good to have an explicit higher-level
version of the C implementation, for example in Maxima, that one could
test against and understand more easily.

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

regarding error estimates in the special functions:

On Mon, 2009-09-07 at 13:13 +0100, Brian Gough wrote:
> 
> I have found the error estimates to be pretty useful when debugging,
> and they are only a small extra cost.

It's not about cycles. Nobody is counting cycles (yet).

It's about the intellectual cost of having all that junk in there. 
It takes up valuable real estate, makes it very difficult for
third party contributors to work with, and it is not based on
any rational methodology.

There is a long discussion about this on the mailing list,
from some years ago. It arose again recently. See the
archived thread starting with:
  http://sourceware.org/ml/gsl-discuss/2004-q3/msg00059.html


Now about cycles. Whether or not the cycle count shows a real
performance issue on current hardware is an open question. I'm
sure it depends on the details, and probably on hardware in
some cases.

But whether or not it is real, there is perceived problem.
See the comments in the above thread. Perceived problems
are real problems, when you are trying to expand the user base.


> Some of the hypergeometric functions are difficult to understand in
> places.  

I wrote them, and I find them nearly impossible to understand.
And I never use them; I would rather spend time finding a more
reliable solution in any specific circumstance. Eventually they
will be scrapped, or at least relegated to some closet in the
library where we store messy things.

As I have said before, having those functions living right
next door to gold-plated Bessel function implementations
is incongruous. This is a _real_ problem because it clouds
the issue of reliability for the whole module.


> Here it would be good to have an explicit higher-level
> version of the C implementation, for example in Maxima, that one could
> test against and understand more easily.

Maybe. If somebody wanted to do that. But it's never that simple.
That would only address some fraction of the possible problems;
the details would remain too different. A combination of methods
is needed.

The first real step would be to refactor so that all the commonality
in methods is better exposed. There are too many places where the
same continued fraction method or series summation method is
being applied, with only slight variation. But I digress.

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Tue, 15 Sep 2009 18:47:01 -0600,
Gerard Jungman wrote:
> The first real step would be to refactor so that all the commonality
> in methods is better exposed. There are too many places where the
> same continued fraction method or series summation method is
> being applied, with only slight variation. But I digress.

If that can be done without changes to the external interface them I
am certainly in favour of it.

-- 
Brian 

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Wed, 9 Sep 2009 16:18:54 +0100,
Jonathan Underwood wrote:
> It'd be nice if it was possible to implement an approach such as
> CESTAC/CADNA in GSL ..... however, I don't see how to do it in C, as
> it seems to require overloading of the basic operators (+/-* etc).

Yes that would be the ideal way to do it, with the error estimate
determined automatically, but as you suggest C is not capable of
handling that--it would have to be done by code generation or a parser
that preprocesses the source (which is how this sort of thing used to
be done in FORTRAN 77 I believe).  

I found the SIAM book "Accuracy and reliability in scientific
computing" a useful reference on this subject, although it is a few
years old now so it may not cover the latest developments.

-- 
Brian Gough

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

>    Linear algebra is a foundational tool in numerical computing.
>    Users need it; other modules in GSL need it. We need a serious
>    solution for this mess.

As I mentioned in previous emails, there isn't a feasible general
solution for LAPACK with the column-major matrix restriction, due to
the combinatorial explosion of cases that must be handled.

We could adopt FLAME, which is a much more general framework, C-based
and faster than LAPACK.  I really think this is a better way to go
than LAPACK.  Unfortunately it doesn't have so many routines at the
moment.

-- 
Brian Gough

Re: GSL 2.0 roadmap (one man's view)

[Rhys Ulerich]

> We could adopt FLAME, which is a much more general framework, C-based
> and faster than LAPACK.  I really think this is a better way to go
> than LAPACK.  Unfortunately it doesn't have so many routines at the
> moment.

+1 on the idea.  The FLAME developers would probably be willing to
fill in the gaps if it meant having FLAME underneath GSL 2.0.  They're
nice folks, and they love displacing LAPACK.

- Rhys

Re: GSL 2.0 roadmap (one man's view)

[Robert G. Brown]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 1003 bytes --]

On Mon, 7 Sep 2009, Rhys Ulerich wrote:

>> We could adopt FLAME, which is a much more general framework, C-based
>> and faster than LAPACK.  I really think this is a better way to go
>> than LAPACK.  Unfortunately it doesn't have so many routines at the
>> moment.
>
> +1 on the idea.  The FLAME developers would probably be willing to
> fill in the gaps if it meant having FLAME underneath GSL 2.0.  They're
> nice folks, and they love displacing LAPACK.

I've never been that fond of LAPACK anyway -- it's not like it is easy
to use and perfectly intuitive anyway.

The biggest issue would/will probably be rationalizing the views of
vector and matrix so they are sufficiently portable and easy to e.g.
pass in and out of ODE solvers and everything else consistently.

    rgb

>
> - Rhys
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Mon, 2009-09-07 at 16:09 +0100, Brian Gough wrote:
> 
> As I mentioned in previous emails, there isn't a feasible general
> solution for LAPACK with the column-major matrix restriction, due to
> the combinatorial explosion of cases that must be handled.


I don't understand "combinatorial explosion". LAPACK only handles
column-major, so it can only be fed column-major. If you want to
use it, all your data structures are born column-major, or you
convert them as needed.

We need to simply allow for the existence of the column-major
world. It's only fair. And it's a simple variation to encode
in any of our data structures.

Every serious C++ vector/matrix library allows for both column-major
and row-major layout, because they know people need it. Templates, or
other gymnastics, are not required. You just have to decide that you
want it.


> We could adopt FLAME, which is a much more general framework, C-based
> and faster than LAPACK.  I really think this is a better way to go
> than LAPACK.  Unfortunately it doesn't have so many routines at the
> moment.

For me, the string "LAPACK" stands for any existing code base that
people will want/need to use. FLAME is just another "LAPACK" for these
purposes. There is nothing for us to "adopt". It's the end user who
"adopts" a package to use. We just have to make it easy for them to
adopt whatever they like.

And if FLAME, or any other library, uses data structures which make it
hard for people to talk to or translate to, then they have made a
mistake.


On Mon, 2009-09-07 at 14:21 -0400, Robert G. Brown wrote:
> I've never been that fond of LAPACK anyway -- it's not like it is easy
> to use and perfectly intuitive anyway.

I know. But it works. And people want it. It's crazy to ignore it.

--
G. Jungman

column-major

[Z F]

--- On Tue, 9/15/09, Gerard Jungman <jungman@lanl.gov> wrote:

> We need to simply allow for the existence of the
> column-major
> world. It's only fair. And it's a simple variation to
> encode
> in any of our data structures.
> 
> Every serious C++ vector/matrix library allows for both
> column-major
> and row-major layout, because they know people need it.
> Templates, or
> other gymnastics, are not required. You just have to decide
> that you
> want it.
> 
>

Could someone, please, explain to me the difference between
column and row major data storage.

I understand it that a multi-dim data is stored in a memory block
which is accessed via

A[i,j,k] = (k*Nj+j)*Ni+i
A[i,j,k] = (i*Nj+j)*Nk+k

Is this correct? or there is more to it? which is which?
If this is correct, what is the big deal with this?
Package which supports one, supports the other automatically.
Left multiplication becomes right multiplication and vice versa. Am I wrong?

Thanks

ZF


      

Re: column-major

[Brian Gough]

At Thu, 17 Sep 2009 20:51:30 -0700 (PDT),
Z F wrote:
> Package which supports one, supports the other automatically.
> Left multiplication becomes right multiplication and vice versa. Am I wrong?

There is a good introduction to the problem in the document "Legacy
BLAS: C Interface to the Legacy BLAS" (section B2.12) at
http://www.netlib.org/blas/blast-forum/

-- 
Brian Gough

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Thu, 27 Aug 2009 17:15:39 -0600,
Gerard Jungman wrote:
>    The GSL native fft implementations suffer in the same way that
>    the GSL native linear algebra implementations suffer. Better
>    solutions are available, from almost all points of view, and
>    the GSL implementations therefore detract from the GSL
>    cachet as a whole. "Serious users should use FFTW" is not
>    a good tagline for the sublibrary; see the discussion of
>    linear algebra above.

In practice it seems like any solution is going come down to saying
people should use FFTW (or LAPACK) one way or another.  Either we
could bundle it, or -- if we don't supply the existing GSL routines --
they could be forced to install it.  Either way it doesn't seem like
the result is actually much different from now.

While getting packages via RPM/DEB is typically fine for desktop use,
I think the ease of doing this in the wider world generally is
underestimated:

- if someone needs a version not available as an RPM/DEB
- on other operating systems
- for cross-compilation

Consider an example of targetting an embedded system for
data-collection with an ARM cpu where the vendor provides only an
older version of GCC, without fortran.  Based on my experience, this
is a realistic scenario in industry.

Currently GSL could be used easily in this scenario via
cross-compilation.  If we depended on external packages it would be
much more difficult.

Re: GSL 2.0 roadmap (one man's view)

[Gerard Jungman]

On Mon, 2009-09-07 at 14:57 +0100, Brian Gough wrote:
> 
> In practice it seems like any solution is going come down to saying
> people should use FFTW (or LAPACK) one way or another.  Either we
> could bundle it, or -- if we don't supply the existing GSL routines --
> they could be forced to install it.  Either way it doesn't seem like
> the result is actually much different from now.

What we type may not be much different. What we say could be quite
different. We need to catalog the use cases for these things. It
mainly comes down to inter-operability of data structures.

In the end, everybody's data structures amount to something plastered
over the top of double *, with some meta-information for layout and
life-cycle control. This includes things likes the row-major vs
column-major distinction, which we should not ignore any longer. 

But this rapidly turns into a discussion of gsl_vector and friends,
which we should save for later. The distinct point is that the
GSL stop-gap implementations are hindering us. I have seen
benchmark web pages, regarding ffts or linear algebra,
where GSL looks like a turd. Well... if the dunce cap fits...


> While getting packages via RPM/DEB is typically fine for desktop use,
> I think the ease of doing this in the wider world generally is
> underestimated:

Packages are the best thing to happen in computer
organization and configuration ever. Period. But anyway...


> Consider an example of targetting an embedded system for
> data-collection with an ARM cpu where the vendor provides only an
> older version of GCC, without fortran.  Based on my experience, this
> is a realistic scenario in industry.

We could have a poll. How many people reading this mailing list
are targeting a platform like this, and how many are sitting
in front of a linux box right now, trying to get their crap
to work, so that they can get on with life?

The fft and linear algebra code need not be sent down the
garbage chute. It will be there for fortran-less ARM prisoners.
But the rest of us should get something better by default.

--
G. Jungman

Re: GSL 2.0 roadmap (one man's view)

[Tuomo Keskitalo]

On 09/16/2009 03:48 AM, Gerard Jungman wrote:

> In the end, everybody's data structures amount to something plastered
> over the top of double *, with some meta-information for layout and
> life-cycle control. This includes things likes the row-major vs
> column-major distinction, which we should not ignore any longer. 
> 
> But this rapidly turns into a discussion of gsl_vector and friends,
> which we should save for later.

Data structures are sort of foundational, shouldn't we discuss them also 
now?

I'd like to know how the developers of GSL ended up using gsl_block 
within gsl_vector and gsl_matrix. Why is it needed?

tda in gsl_matrix struct is a bit of a mystery to me, what is it used for?

I agree that use of other libraries and even other languages with GSL 
should be made as easy as possible. Data structures with dimension > 1 
should be explicit about the data being stored as "column-major" or 
"row-major".

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: GSL 2.0 roadmap (one man's view)

[Robert G. Brown]

On Sun, 20 Sep 2009, Tuomo Keskitalo wrote:

> I agree that use of other libraries and even other languages with GSL should 
> be made as easy as possible. Data structures with dimension > 1 should be 
> explicit about the data being stored as "column-major" or "row-major".

I agree completely, and I further strongly suggest that a) the actual
data in any tensor form be stored in a single contiguous vector whose
address is available to the user and not hidden inside an opaque data
object; b) the tensor itself be a **...*pointer, cast to the desired
type, packed with the offsets into this vector; c) that the data should
be stored/packed using the same convention that C already uses to pack
e.g. v[i], m[i][j] declared data types.  Yes, this leaves one with the
perpetual conflict with e.g. fortran libraries, but GSL is a C library,
not a fortran library, and it should conform to the common usages and
standards of C and not try to interpolate to make fortran programmers
comfortable at the expense of C programmers.

My reasons for each suggestion are straightforward.  a) permits many
routines, e.g. ODE solvers, to be written for general purpose tensors by
letting them trivially directly access the vector data associated with
the tensor; b) allows the caller to use meaningful tensor notation to
evaluate functional forms throughout their code instead of using get or
set utilities to access particular elements of the tensor; c) makes the
code portable in critical ways -- it would be silly to have

for(i=0;i<imax;i++)
   for(j=0;j<jmax;j++)
      cmatrix[i][j] = gslmatrix[i][j];

(where cmatrix is a matrix allocated via e.g. double
cmatrix[imax][jmax] and gslmatrix is one created via the process
suggested above) do the wrong thing, and higher rank tensors should
obviously satisfy the same convention to make programming contractions
etc maximally simple.

Since MANY GSL routines could be affected by this -- in particular ODE
solution, linear algebra routines (if/when this is added),
multidimensional quadrature routines (ditto) -- the discussion should
occur at the very beginning.

    rgb

>
> -- 
> Tuomo.Keskitalo@iki.fi
> http://iki.fi/tuomo.keskitalo
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap (one man's view)

[Rhys Ulerich]

> a) the actual data in any tensor form be stored in a single contiguous
> vector whose address is available to the user and not hidden inside an
> opaque data object;
> b) the tensor itself be a **...*pointer, cast to the desired
> type, packed with the offsets into this vector

Given both (a) and (b) everyone pays storage overhead for indexing
twice.  The first time is to store the leading dimensions for (a)-like
access, and the second time is to store the pointers-to-pointers for
(b)-like access.

Since (a) already exists today, if you want (b) it can be provided as
a view.  Your cost will be storage overhead for pointers-to-pointers,
dereferencing delay, and potential cache misses if the
pointers-to-pointers for (b)-like access are far away from (a)'s data.

I think those that want (b) should have a view that accomplishes it,
but that people using (a) only shouldn't pay for the overhead.  I
believe these are the reasons that both the C and Fortran numerics
communities pretty consistently stay towards (a)-like storage under
the covers.

Also, to be fair to the current design, gsl_matrix and gsl_vector are
structs, but they certainly aren't opaque.

- Rhys

Re: GSL 2.0 roadmap (one man's view)

[Robert G. Brown]

On Sun, 20 Sep 2009, Rhys Ulerich wrote:

>> a) the actual data in any tensor form be stored in a single contiguous
>> vector whose address is available to the user and not hidden inside an
>> opaque data object;
>> b) the tensor itself be a **...*pointer, cast to the desired
>> type, packed with the offsets into this vector
>
> Given both (a) and (b) everyone pays storage overhead for indexing
> twice.  The first time is to store the leading dimensions for (a)-like
> access, and the second time is to store the pointers-to-pointers for
> (b)-like access.
>
> Since (a) already exists today, if you want (b) it can be provided as
> a view.  Your cost will be storage overhead for pointers-to-pointers,
> dereferencing delay, and potential cache misses if the
> pointers-to-pointers for (b)-like access are far away from (a)'s data.
>
> I think those that want (b) should have a view that accomplishes it,
> but that people using (a) only shouldn't pay for the overhead.  I

I agree, which is why I suggest that a) always be directly and simply
available for people who are willing to write the somewhat ugly code (or
macros) required to dereference without the additional layer of
overhead and to facilitate other things such as ODE soolutions.

Bear in mind that code performance isn't always the primary issue, even
for numerical programs.  Sometimes ease of programming and/or debugging
are important, or the numerical task takes so little time that nobody
cares about the possibly small increase in execution time.  The
dereferencing code also gets pretty ugly and non-portable for
non-rectangular tensors of high rank with indices that don't always
start at 0.  Even if one's goal is to end up with a fast program with
inline dereferencing an minimal indexing overhead, it might well be VERY
USEFUL to have a direct and straightforward tensor view of the data
vector to use to write an easily debugged version that can in turn be
used to debug a fast version.

> believe these are the reasons that both the C and Fortran numerics
> communities pretty consistently stay towards (a)-like storage under
> the covers.
>
> Also, to be fair to the current design, gsl_matrix and gsl_vector are
> structs, but they certainly aren't opaque.

And they aren't extensible (or at least extended) to tensors of higher
rank.  Scalability matters.

    rgb

>
> - Rhys
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Sun, 20 Sep 2009 10:30:39 -0500,
Rhys Ulerich wrote:
> I think those that want (b) should have a view that accomplishes it,
> but that people using (a) only shouldn't pay for the overhead. 

Yes, it would be reasonable to have a generalised view built on top of
the existing objects like this.

Re: GSL 2.0 roadmap (one man's view)

[Rhys Ulerich]

> tda in gsl_matrix struct is a bit of a mystery to me, what is it used for?

From the manual [1]: "The physical row dimension tda, or trailing
dimension, specifies the size of a row of the matrix as laid out in
memory."

Say A is a (gsl_matrix *) and A->data[offset] is the matrix entry
A_{i,j}, then A->data[offset+A->tda] is the matrix entry A_{i+1,j}.
Note that both tda and size2 are necessary because our matrix A may be
a view of some larger matrix.

- Rhys

[1] http://www.gnu.org/software/gsl/manual/html_node/Matrices.html

Re: GSL 2.0 roadmap (one man's view)

[Brian Gough]

At Sun, 20 Sep 2009 12:35:59 +0300,
Tuomo Keskitalo wrote:
> I'd like to know how the developers of GSL ended up using gsl_block 
> within gsl_vector and gsl_matrix. Why is it needed?

The block type corresponds to the C++ valarray class (the gsl_block
struct members are the same as those in C++ valarrays).  The idea was
that it would simplify interoperability with C++.

-- 
Brian Gough

Re: GSL 2.0 roadmap

[Brian Gough]

At Fri, 21 Aug 2009 21:41:53 +0100,
Brian Gough wrote:
> I plan to make a branch for release 2.0 after the next release 1.13
> which will hopefully be next week. 

I've made the 2.0 tarball.  Just sorting some problems with support
for multiple GPG signatures on ftp.gnu.org before it can go out.  I'm
away next week so I'll probably upload it around 7 September.