Matvec.module

(* 

The matrix and vector classes were originally developed by S.
Omohundro.  Later, however, Matt Kennel put a tremendous effort into
defining new abstractions and links into the BLAS libraries.  These
are his notes on the fundamental numerical matrix classes for Sather.

First issue to clarify: what is the difference between a "matrix" and
a "two-dimensional array"?

A: A matrix is a finite-dimensional linear operator over a given
coordinate system (the vector space).  As such it pretty much only
makes sense to have matrices whose elements are ring elements, that
is, can add, subtract and multiply, have additive and and
multiplicative identities, as well as additive inverses.

 At this stage, I'm going to require in addition that the matrix
elements also be _field elements_, requiring that they have
multiplicative inverses and hence division.  The only fields that I
can really see right now being relevant are the real and complex
numbers, each implemented as single precision or double precision
floating point.  Maybe quaternions would be useful for somebody, but
one must keep in mind that quaternions are not Abelian (commutative)
and so care would have to be taken to define proper order of
operations.

 Matrices are often _implemented_ in terms of two-dimensional arrays
of numbers, but this is not fundamental.  Indeed in realistic
numerical problems it may make sense to implement matrices in terms of
only that small subset of non-zero elements {sparse matrices} or,
perhaps as a combination of principal components and significant
eigenvectors or whatever.  In linear algebra computations there are
often "factorizations" of matrices: M -> L*U.  A structure consisting
of "L*U" is a good representation for a matrix just as well as a dense
2-d array of numbers.

 Still, the initial goal for the matrix classes is to provide convenient
access to standard computational routines for single and double precision
real and complex matrices.
 
 The Basic matrix classes only include routines for generic dense
rectangular matrices, and no sophisticated linear algebra algorithms.
    
 By contrast the Fancy matrix classes will include linear algebra
computations such as solving systems of linear equations,
eigenproblems, and various kinds of matrix decompositions on dense and
special sparse matrices.

 The Basic matrix class will only assume and require the presence of
BLAS and not LAPACK.  Furthermore, it will NOT provide all of the BLAS
functionality, just those parts relevant to operations on general
dense matrices and vectors.  The Fancy matrix classes will extend the
basic matrix classes by adding most of of the BLAS functionality and
some reasonable subset of LAPACK.  The extent of this depends on how
much gets done by me and other volunteers.
 
 The Basic matrix classes are, in fact, implemented as a
two-dimensional array of numbers.  We choose the *column-major* layout
in memory, following Fortran convention, so as to aid interfacing with
the huge corpus of Fortran numerical routines.  Most other langauges,
e.g. C use row-major arrays by default.  The rather feeble
``built-in'' 2-d arrays in C++ add nothing beyond those in C. C++
matrix classes use a wide variety of representations, but no matrix
class is yet standard.

 The goal of the Sather Basic Matrix Class is to provide a standard
single matrix class design as good as the best C++ classes, with a
standard high-performance implementation as good as the best C++
classes.

 Any successful program library requires a thoughtful design.  It
should provide enough diversity in functionality to support a wide
range of applications, and yet the library should be as simple as
possible.  Better to have a few key mechanisms instead of a profusion
of similar, but not identical mechanisms.  Finally, it should provide
maximum performance.
 
 The tension between these three is eternal and fundamental.  Its
resolution requires creativity and experience.
 
 Fortunately, the computational world has extensive experience with
linear algebra, and thus the difficult intellectual classifications
have already been done.

 In particular, the Basic Linear Algebra Subroutines (BLAS) are the
standard "core" algorithms for contemporary "standard" numerical
subroutines.
 
 In particular, the famous LINPACK and EISPACK, now superseded by
LAPACK, use BLAS routines.
 
 BLAS routines sometimes come written in assembly expertly tuned for
high performance on various platforms.  Even without this, there are
free C and Fortran sources for the BLAS.  {LAPACK and BLAS are
available for download on the Internet from NETLIB, the de facto
standard repository for computational software.  You can use FTP or
the WWW and connect to "netlib.att.com" or "netlib.ornl.gov".}
 
 Thus the standard implementation of Sather's matrix classes will rely
on BLAS routines to for the computational gruntwork.  The fuller
matrix classes will use BLAS and LAPACK to provide easy-to-use
interfaces to a variety of sophisticated high-performance numerical
algorithms.
 
 Out of the three factors mentioned above, the BLAS routines have been
developed mostly for performance.  Their algorithmic interfaces to the
human programmer are not particularly friendly.  The Sather Basic
matrix classes are intended to provide a convenient front-end with
idiomatic Sather functionality to traditional Fortran routines without
substantially sacrificing performance.  In choosing member routines, I
put a higher priority on providing many explicit variations of
operations for the sake of having the least overhead instead of
accepting some run-time penalty in favor of having a tighter and
smaller class interface.

 This is similar to the goal of the LAPACK++ library for C++.  I think
this is the best C++ matrix class library, and I will keep it in mind
while designing the Sather library, though I do not intend to copy it
slavishly.  In particular LAPACK++ pretty much requires LAPACK as well
as BLAS; I hope to have a basic subset which does not LAPACK to
compile.  The LAPACK++ library accesses individual elements through
double indirection; its matrix classes contain a pointer to
dynamically allocated raw storage.  That allows 'matrix references'
that denote subsets of existing matrices (3rd through 5th row, 2nd
through 9th column of A) to be the same class as the matrix class
itself simplifying the interface at some performance cost.  I choose
not to pursue this design for the sake of performance.  As is natural
for Sather, the heap storage for the matrix will be contained "in" the
object itself via an inheritance path to AREF{T}.  At some point
'matrix reference' classes may be written if there is a need.
 
 The Sather compiler itself has optimization technology capable of
making naive executable code for numerical kernels that is much faster
than naive C++ and C, though handtuned C can achieve somewhat better
results.  
 
 However, for today's microprocessors, kernels hand tuned for a
particular architecture and cache can perform far better than even
fairly careful, portable C. Most users will be better served by front
ends to those fast and reliable algorithms in BLAS and LAPACK rather
than new development in Sather.
 
 Other sorts of algorithms, e.g. minimizations, ODE integrations,
genetic algorithms and the like would benefit more from the modern
data structuring and object-orientation facilities of Sather, and hence
native rewrites would be more rewarding.  I hope somebody tries to
work on these, as I don't have the time or expertise.

*)

abs_mat.sa  -has abs_mat.sa $MAT
mat.sa  -has mat.sa MAT NUMERIC_MAT
matd.sa -has matd.sa MATD
matcpx.sa -has matcpx.sa MATCPX MATCPXD
abs_vec.sa -has abs_vec.sa $VEC
vec.sa  -has vec.sa VEC VEC_LENGTH_MIXIN 
vecd.sa -has vecd.sa VECD
veccpx.sa -has veccpx.sa VECCPX_LENGTH_MIXIN VECCPX VECCPXD
test_matvec.sa -has test_matvec.sa TEST_MATVEC