gsl-design.texi

开放gsl矩阵运算

unmaintainable code.  However, it's easy to go down that road if you
don't think about it ahead of time.

It is better to choose simplicity over completeness.  In designing new
parts of the library keep modules independent where possible. If
interdependencies between modules are introduced be sure about where you
are going to draw the line.

@node Code Reuse, Standards and conventions, Design of  Numerical Libraries, Design
@section Code Reuse

It is useful if people can grab a single source file and include it in
their own programs without needing the whole library.  Try to allow
standalone files like this whenever it is reasonable.  Obviously the
user might need to define a few macros, such as GSL_ERROR, to compile
the file but that is ok.  Examples where this can be done: grabbing a
single random number generator.

@node Standards and conventions, Background and Preparation, Code Reuse, Design
@section Standards and conventions

The people who kick off this project should set the coding standards and
conventions.  In order of precedence the  standards that we follow are,

@itemize @bullet
@item We follow the GNU Coding Standards.
@item We follow the conventions of the ANSI Standard C Library.
@item We follow the conventions of the GNU C Library.
@item We follow the conventions of the glib GTK support Library.
@end itemize

The references for these standards are the @cite{GNU Coding Standards}
document, Harbison and Steele @cite{C: A Reference Manual}, the
@cite{GNU C Library Manual} (version 2), and the Glib source code.

If the project has a philosophy it is to "Think in C".  Since we are
working in C we should only do what is natural in C, rather than trying
to simulate features of other languages.  If there is something which is
unnatural in C and has to be simulated then we avoid using it. If this
means leaving something out of the library, or only offering a limited
version then so be it.  It is not worthwhile making the library
over-complicated.  There are numerical libraries in other languages, and
if people need the features of those languages it would be sensible for
them to use the corresponding libraries, rather than coercing a C
library into doing that job.  

It should be borne in mind at all time that C is a macro-assembler.  If
you are in doubt about something being too complicated ask yourself the
question "Would I try to write this in macro-assembler?" If the answer
is obviously "No" then do not try to include it in GSL. [BJG]

It will be useful to read the following paper,

@itemize @asis
@item 
Kiem-Phong Vo, ``The Discipline and Method Architecture for Reusable
Libraries'', Software - Practice & Experience, v.30, pp.107-128, 2000.
@end itemize
@noindent
It is available from
@url{http://www.research.att.com/sw/tools/sfio/dm-spe.ps} or the earlier
technical report Kiem-Phong Vo, "An Architecture for Reusable Libraries"
@url{http://citeseer.nj.nec.com/48973.html}.

There are associated papers on Vmalloc, SFIO, and CDT which are also
relevant to the design of portable C libraries.

@itemize @asis
@item
Kiem-Phong Vo, ``Vmalloc: A General and Efficient Memory
Allocator''. Software Practice & Experience, 26:1--18, 1996.

@url{http://www.research.att.com/sw/tools/vmalloc/vmalloc.ps}
@item
Kiem-Phong Vo. ``Cdt: A Container Data Type Library''. Soft. Prac. &
Exp., 27:1177--1197, 1997

@url{http://www.research.att.com/sw/tools/cdt/cdt.ps}
@item
David G. Korn and Kiem-Phong Vo, ``Sfio: Safe/Fast String/File IO'',
Proceedings of the Summer '91 Usenix Conference, pp.  235-256, 1991.

@url{http://citeseer.nj.nec.com/korn91sfio.html}
@end itemize


@node Background and Preparation, Documentation, Standards and conventions, Design
@section Background and Preparation

Before implementing something be sure to research the subject
thoroughly.  This will save a lot of time in the long-run.  The two most
important steps are,

@enumerate
@item
to determine whether there is already a free library (GPL or
GPL-compatible) which does the job.  If so, there is no need to
reimplement it.  Carry out a search on Netlib, GAMs, na-net,
sci.math.num-analysis and the web in general. This should also provide a
list of existing proprietary libraries which are relevant.

@item 
make a comparative survey of existing implementations in the
commercial/free libraries. Examine the typical APIs, methods of
communication between program and subroutine, and classify them so that
you are familiar with the key concepts or features that an
implementation may or may not have, depending on the relevant tradeoffs
chosen.  Be sure to review the documentation of existing libraries for
useful references.

@item
read up on the subject and determine the state-of-the-art.  Find the
latest review papers.  A search of the following journals should be
undertaken.

@itemize @asis
@item ACM Transactions on Mathematical Software
@item Numerische Mathematik
@item Journal of Computation and Applied Mathematics
@item Computer Physics Communications
@item SIAM Journal of Numerical Analysis
@item SIAM Journal of Scientific Computing
@end itemize
@end enumerate

@node Documentation, Namespace, Background and Preparation, Design
@section Documentation
Documentation: the project leaders should give examples of how things
are to be documented.  High quality documentation is absolutely
mandatory, so documentation should introduce the topic, and give careful
reference for the provided functions. The priority is to provide
reference documentation for each function. It is not necessary to
provide tutorial documentation.

Use free software, such as GNU Plotutils, to produce the graphs in the
manual.

Some of the graphs have been made with gnuplot which is not truly free
(or GNU) software, and some have been made with proprietary
programs. These should be replaced with output from GNU plotutils.

@node Namespace, Header files, Documentation, Design
@section Namespace

Use @code{gsl_} as a prefix for all exported functions and variables.

Use @code{GSL_} as a prefix for all exported macros.

All exported header files should have a filename with the prefix @code{gsl_}.

All installed libraries should have a name like libgslhistogram.a

Any installed executables (utility programs etc) should have the prefix
@code{gsl-} (with a hyphen, not an underscore).

@node Header files, Target system, Namespace, Design
@section Header files

Installed header files should be idempotent, i.e. surround them by the
preprocessor conditionals like the following,

@example
#ifndef __GSL_HISTOGRAM_H__
#define __GSL_HISTOGRAM_H__
...
#endif /* __GSL_HISTOGRAM_H__ */
@end example

@node Target system, Function Names, Header files, Design
@section Target system

The target system is ANSI C, with a full Standard C Library, and IEEE
arithmetic.

@node Function Names, Object-orientation, Target system, Design
@section Function Names

Each module has a name, which prefixes any function names in that
module, e.g. the module gsl_fft has function names like
gsl_fft_init. The modules correspond to subdirectories of the library
source tree.

@node Object-orientation, Comments, Function Names, Design
@section Object-orientation

The algorithms should be object oriented, but only to the extent that is
easy in portable ANSI C.  The use of casting or other tricks to simulate
inheritance is not desirable, and the user should not have to be aware
of anything like that.  This means many types of patterns are ruled
out.  However, this is not considered a problem -- they are too
complicated for the library.  

Note: it is possible to define an abstract base class easily in C, using
function pointers.  See the rng directory for an example. 

When reimplementing public domain fortran code, please try to introduce
the appropriate object concepts as structs, rather than translating the
code literally in terms of arrays.  The structs can be useful just
within the file, you don't need to export them to the user.

For example, if a fortran program repeatedly uses a subroutine like,

@example
SUBROUTINE  RESIZE (X, K, ND, K1)
@end example
@noindent
where X(K,D) represents a grid to be resized to X(K1,D) you can make
this more readable by introducing a struct,

@example
struct grid @{
    int nd;  /* number of dimensions */
    int k;   /* number of bins */
    double * x;   /* partition of axes, array of size x[k][nd] */
@}

void
resize_grid (struct grid * g, int k_new)
@{
...
@}
@end example


@node Comments, Minimal structs, Object-orientation, Design
@section Comments

Follow the GNU Coding Standards.  A relevant quote is,

``Please write complete sentences and capitalize the first word.  If a
lower-case identifier comes at the beginning of a sentence, don't
capitalize it!  Changing the spelling makes it a different identifier.
If you don't like starting a sentence with a lower case letter, write
the sentence differently (e.g., "The identifier lower-case is ...").''

@node Minimal structs, Algorithm decomposition, Comments, Design
@section Minimal structs

We prefer to make structs which are minimal.  For example, if a certain
type of problem can be solved by several classes of algorithm (e.g. with
and without derivative information) it is better to make separate types
of struct to handle those cases.  i.e. run time type identification is
not desirable.

@node Algorithm decomposition, Memory allocation and ownership, Minimal structs, Design
@section Algorithm decomposition

Iterative algorithms should be decomposed into an INITIALIZE, ITERATE,
TEST form, so that the user can control the progress of the iteration
and print out intermediate results.  This is better than using
call-backs or using flags to control whether the function prints out
intermediate results.  In fact, call-backs should not be used -- if they
seem necessary then it's a sign that the algorithm should be broken down
further into individual components so that the user has complete control
over them.  

For example, when solving a differential equation the user needs to be
able to advance the solution by individual steps, while tracking a
realtime process.  This is only possible if the algorithm is broken down
into step-level components.  Higher level decompositions would not give
sufficient flexibility.

@node Memory allocation and ownership, Memory layout, Algorithm decomposition, Design
@section Memory allocation and ownership

Functions which allocate memory on the heap should end in _alloc
(e.g. gsl_foo_alloc) and be deallocated by a corresponding _free function
(gsl_foo_free).

Be sure to free any memory allocated by your function if you have to
return an error in a partially initialized object.

Don't allocate memory 'temporarily' inside a function and then free it
before the function returns.  This prevents the user from controlling
memory allocation.  All memory should be allocated and freed through
separate functions and passed around as a "workspace" argument.  This
allows memory allocation to be factored out of tight loops.

@node Memory layout, Linear Algebra Levels, Memory allocation and ownership, Design
@section Memory layout

We use flat blocks of memory to store matrices and vectors, not C-style
pointer-to-pointer arrays.  The matrices are stored in row-major order
-- i.e. the column index (second index) moves continuously through memory. 

@node Linear Algebra Levels, Exceptions and Error handling, Memory layout, Design
@section Linear Algebra Levels

Functions using linear algebra are divided into two levels:

For purely "1d" functions we use the C-style arguments (double *,
stride, size) so that it is simpler to use the functions in a normal C
program, without needing to invoke all the gsl_vector machinery.

The philosophy here is to minimize the learning curve. If someone only
needs to use one function, like an fft, they can do so without having
to learn about gsl_vector.  

This leads to the question of why we don't do the same for matrices.
In that case the argument list gets too long and confusing, with
(size1, size2, tda) for each matrix and potential ambiguities over row
vs column ordering. In this case, it makes sense to use gsl_vector and
gsl_matrix, which take care of this for the user.

So really the library has two levels -- a lower level based on C types
for 1d operations, and a higher level based on gsl_matrix and
gsl_vector for general linear algebra.

Of course, it would be possible to define a vector version of the
lower level functions too. So far we have not done that because it was
not essential -- it could be done but it is easy enough to get by
using the C arguments, by typing v->data, v->stride, v->size instead.
A gsl_vector version of low-level functions would mainly be a
convenience.

@node Exceptions and Error handling, Persistence, Linear Algebra Levels, Design
@section Exceptions and Error handling

The basic error handling procedure is the return code (see gsl_errno.h
for a list of allowed values).  Use the GSL_ERROR macro to mark an
error.  The current definition of this macro is not ideal but it can be
changed at compile time. 

You should always use the GSL_ERROR macro to indicate an error, rather
than just returning an error code. The macro allows the user to trap
errors using the debugger (by setting a breakpoint on the function
gsl_error).  

The only circumstances where GSL_ERROR should not be used are where the
return value is "indicative" rather than an error -- for example, the
iterative routines use the return code to indicate the success or
failure of an iteration. By the nature of an iterative algorithm
"failure" (a return code of GSL_CONTINUE) is a normal occurrence and
there is no need to use GSL_ERROR there.

Be sure to free any memory allocated by your function if you return an
error (in particular for errors in partially initialized objects).

@node Persistence, Using Return Values, Exceptions and Error handling, Design
@section Persistence

If you make an object foo which uses blocks of memory (e.g. vector,
matrix, histogram) you can provide functions for reading and writing
those blocks,

@example
int gsl_foo_fread (FILE * stream, gsl_foo * v);
int gsl_foo_fwrite (FILE * stream, const gsl_foo * v);
int gsl_foo_fscanf (FILE * stream, gsl_foo * v);
int gsl_foo_fprintf (FILE * stream, const gsl_foo * v, const char *format);
@end example
@noindent
Only dump out the blocks of memory, not any associated parameters such
as lengths.  The idea is for the user to build higher level input/output
facilities using the functions the library provides.  Use the functions

@example
int gsl_block_fread (FILE * stream, gsl_block * b);
int gsl_block_fwrite (FILE * stream, const gsl_block * b);
int gsl_block_fscanf (FILE * stream, gsl_block * b);
int gsl_block_fprintf (FILE * stream, const gsl_block * b, const char *format);
@end example
@noindent
or

@example
int gsl_block_raw_fread (FILE * stream, double * b, size_t n, size_t stride);
int gsl_block_raw_fwrite (FILE * stream, const double * b, size_t n, size_t stri
de);
int gsl_block_raw_fscanf (FILE * stream, double * b, size_t n, size_t stride);
int gsl_block_raw_fprintf (FILE * stream, const double * b, size_t n, size_t str
ide, const char *format);
@end example
@noindent
to do the actual reading and writing.