faq_meschach.txt

C语言版本的矩阵库

		Meschach Frequently Asked Questions
			    Version 0.1
			  2nd March 1994


This is a list of frequently asked questions regarding Meschach.
It is maintained by David Stewart and Zbigniew Leyk, the developers and
maintainers of Meschach.  A copy of this is kept in the anonymous ftp
area

    thrain.anu.edu.au : /pub/meschach/FAQ.meschach

----------------------------------------------------------------------

Q1.  What is Meschach?
Q2.  Why not use C versions of LINPACK, EISPACK or LAPACK?
Q2a. What are the differences of Meschach from LAPACK or EISPACK in C?
Q2b. How reliable are the routines compared to the above packages?
Q3.  How do I get documentation?
Q4.  I've found a bug in Meschach!  What do I do?
Q5.  It won't compile on my BrandX machine running X-OS.  What do I do?
Q6.  Is there a C++ version of Meschach?
Q7.  Meschach doesn't have an xxxx routine.  Why not?  What do I do?
Q8.  Can I use Meschach with MATLAB?
Q9.  I don't want ALL of Meschach.  How do I get just a subset?
Q10. What does Meschach mean?
Q11. My code that uses Meschach seems to run very slowly/uses incredible
	amounts of memory.  Why?
Q12. I have a huge problem and I keep running out of memory.  What do I do?
Q13. How did you come to write Meschach?
Q14. Complex part of Meschach
Q15. Can I do an index search in order to know what routines exist?

----------------------------------------------------------------------

A1.  What is Meschach?

Meschach is a library of functions written in C for performing matrix
computations.  The operations can be (briefly) summarised as

	* linear combinations of vectors and matrices; inner products;
	  matrix-vector and matrix-matrix products
	* solving linear systems of equations
	* solving least-squares problems
	* computing eigenvalues (or singular values)
	* finding "nice" bases for subspaces like ranges of matrices and
	  null spaces.

There are, however, a lot of different ways of, for example, solving a
system of linear equations depending on whether the matrix is real or
complex, symmetric, positive definite or indefinite, sparse or has some
other structure.  The rule in computational mathematics is: Exploit
Structure.  So there is the standard Gaussian Elimination method for
general real (and general complex) matrices, as well as Cholesky
factorisation (for positive definite symmetric matrices), which is about
twice as fast, LDL factorisation which just needs symmetric matrices, but
is not always well-behaved, BKP (Bunch-Kaufman-Parlett) factorisation for
indefinite symmetric matrices which is well-behaved but more complex,
QR factorisation which is best used for least squares problems.
And then there are sparse matrices, and a variety of different techniques
for different types of sparse matrices.

A list of the features of Meschach follows.

	o  dynamic allocation, de-allocation, re-sizing and copying of objects
	o  dense complex matrices and vectors as well as real matrices
		and vectors
	o  input and output routines for these objects, and MATLAB
		save/load format
	o  error/exception handling
	o  basic numerical linear algebra -- linear combinations, inner
		products, matrix-vector and matrix-matrix products
		including transposed and adjoint forms
	o  vector min, max, sorting, componentwise products, quotients
	o  dense matrix factorise and solve -- LU, Cholesky, LDL^T, QR,
		QR with	column pivoting, symmetric indefinite (BKP)
	o  dense matrix factorisation update routines -- LDL^T, QR
		(real matrix updates only)
	o  eigenvector/eigenvalue routines -- symmetric, real Schur
		decomposition, SVD, extract eigenvector
	o  sparse matrix "utility" routines
	o  sparse matrix factorise and solve -- Cholesky, LU and BKP
		(Bunch-Kaufman-Parlett symmetric indefinite factorisation)
	o  sparse incomplete factorisations -- Cholesky and LU
	o  iterative techniques -- pre-conditioned conjugate gradients,
		CGNE, LSQR, CGS, GMRES, MGCR, Lanczos, Arnoldi
	o  allowance for "procedurally defined" matrices in the iterative
		techniques
	o  various "torture" routines for checking aspects of Meschach
	o  memory tracking for locating memory leaks

Here is a sample piece of code using Meschach to read in a matrix A and a
vector b and to solve A.x = b for x and print the result.

MAT  *A, *Acopy;
VEC  *b, *x;
PERM *pivot;

A = m_input(MNULL);    /* read in matrix from stdin */
Acopy = m_copy(A,MNULL);  /* create & copy matrix */
b = v_input(VNULL);    /* read in vector from stdin */
pivot = px_get(A->m);  /* get pivot permutation -- size = # rows of A */
LUfactor(A,pivot);     /* LU factorisation */
x = LUsolve(A,pivot,b,VNULL);  /* ... and solve A.x = b; result is in x */
/* .... and print out solution with a comment */
printf("# x =\n");    v_output(x);

Meschach aims to be an easy-to-use, easy-to-debug-your-code library.
The basic data objects (vectors, matrices etc.) are represented by
self-contained data structures which can be created, destroyed and resized
at will.  Workspace arrays do not need to be passed to routines (except
for some internal routines).  Operations are continually checked to see if
they are valid so that errors are quickly located.  Unlike matrix
interpreter/calculators such as MATLAB, the user is given control over
memory usage and can efficiently re-use memory which requires garbage
collection in a more interpretive environment.
----------------------------------------------------------------------

A2.  Why not use C versions of LINPACK, EISPACK or LAPACK?

LINPACK, EISPACK and LAPACK are well-known Fortran-language packages of
routines for a variety of problems in matrix computations.  With the magic
of translators like f2c, why not use the translated versions of these
libraries?  Well, a number of people are already doing this; after all,
LINPACK etc. are public domain and they have been incorporated in public
domain software such as Octave and RlaB, to name but two.

But, as with all translations from Fortran, these do not use any of the
facilities that Fortran lacks.  In particular, there is no attempt to
structure the data, nor is there any use of dynamic memory allocation,
C's error handling (using setjmp() and longjmp()).  Finally there is no
input/output.

Meschach uses not only the data structuring facilities in C to provide
self-contained data structures, but this is used with the error handling to
ensure, for example, that the matrices passed to a matrix multiplication
routine have compatible sizes.  The data structuring is also used with the
dynamic memory allocation to allow efficient resizing -- vectors are only
physically reallocated if the requested size is above the current "high
water mark".  Input and output are both handled using C's I/O facilities to
provide output that is both machine and human-readable.  Binary
MATLAB-compatible output (for dense matrices) is also done.  And another
useful aspect of Meschach's I/O is that it allows for comments in the data.

As well, none of the above packages has anything to do with sparse
matrices.  There are, of course, sparse matrix packages that are freely
available, but even there, they do not seem to handle, in an integrated
way, both sparse and dense matrices.  Meschach has not only sparse matrix
data structures and the routines to make use of them, but as the sparse
matrices are self-contained data structures, use of sparse and dense
matrices can be intermixed.

You may also like to look at answer A2a.
----------------------------------------------------------------------

A2a. What is the difference of meschach from LAPACK or EISPACK in C?

First: The biggest difference is that Meschach is written from scratch in
C, and works in terms of data structures that are self-contained
mathematical objects (matrices, vectors, permutations, sparse matrices
etc).  These can be created, destroyed and re-sized at will.

Second: LAPACK and EISPACK have no sparse matrix routines, and no iterative
routines; Meschach has both and includes sparse factorisations, both with
fill-in and "incomplete" factorisations.

Third: The algorithms are similar; we have **NOT** taken
LINPACK/EISPACK/LAPACK and re-implemented this in C.

Fourth: Meschach makes comprehensive use of C's facilities (as well as
structures a.k.a. records) such as error handling, and dynamic memory
allocation.  E.g. workspace arrays do not need to be passed.  They are
created and re-sized as necessary where they are used.  With Meschach 1.2
you can destroy any of these local workspace vectors etc if memory is
short.

Fifth: Meschach is smaller in the amount of memory needed, both for source
and executable than e.g. LAPACK.  Full source code in shar files is under
1Mbyte for Meschach, more like 4-5Mbyte for LAPACK.  There **are** a number
of things in LAPACK that are not in Meschach (e.g. banded complex
factorisations), but then again, there is a package of sparse matrix
operations in Meschach that is not in LAPACK.  We suspect that it is partly
due to the way the software has been constructed.  Meschach has been a 1 or
2 person effort; LAPACK has had something like 10 full-time programmers.
We've had to be economical in terms of how the software has been developed.
---------------------------------------------------------------------------

A2b. How reliable are the routines compared to the above packages?

Well, LINPACK and EISPACK have been around for nearly 20 years.  The bugs
should be found by now!  So we can't honestly claim that we have yet
achieved **THAT** level of reliability.  But we believe that it probably
approaches that level of reliability for the more commonly used routines.
Because there is a considerable amount of error/exception checking in
Meschach, bugs have been relatively easy to track down.

A number of bugs have been found in Meschach 1.2a (and there will be more).
But these have been mostly minor and easily fixed.
------------------------------------------------------------------------

A3.  How do I get documentation?

There is some basic documentation that comes with the source code of
Meschach and there will be more on that later.
But the comprehensive documentation has been published as a
book (both a users' and a reference manual) by

	Centre for Mathematics and its Applications (CMA)
	School of Mathematical Sciences
	Australian National University
	Canberra, ACT 0200
	Australia

as "Meschach: Matrix Computations in C", volume 32 in the "Proceedings of
the CMA".  The cost is about A$30 (approx. US$22) + postage/handling.
Orders can be passed to the CMA via David Stewart
([e-mail]).  The manual is nearly 250pp long.

If you just want to check out the sort of things that Meschach has before
committing yourself to buying the manual, there is the manual for the
**OLD** version of Meschach which is available over ftp at

    ftpmaths.anu.edu.au : /pub/meschach/version1.1b/bookdvi.tar [.Z or .gz]

You can also have a look at some of the documentation that comes with the
source code.  In the DOC directory there is tutorial.txt which is an ASCII
file containing a tutorial introduction to Meschach.  There is also
fnindex.txt which contains an index (in alphabetical order) of the
functions in Meschach and a brief description of what they do.

The calling sequences can be found from the header files and the actual
source code.  Of course, this is no substitute for a comprehensive manual,
but it does give you a start with using Meschach.  The torture.c and
*tort.c files can also be useful in showing how Meschach can be used.
----------------------------------------------------------------------

A4.  I've found a bug in Meschach!  What do I do?

Congratulations!  Let us know!  Not all bugs are serious --- most of the
bugs found by users are installation problems.
Most bugs in Meschach are very easily fixed.  Please try to determine where
the bug is occuring (which routine, preferably) and what data results in a
problem.  Send that information to us and we will aim to fix the problem as
quickly as possible.  Send the information to either of us.  We will pass
it on to the other.

Sometimes the torture test(s) give a message "Error: ..." indicating that
there may be something wrong.   Please have a look at the size of the error
that is printed out with the "Error:  ..." message.  If this is of the
order of 1e-15 or 1e-14 it is unlikely that there is anything seriously
wrong.  The reason is that torture and the other testing programs generally
use randomly generated problems, which may on occasion be more
ill-conditioned than expected and give slightly larger errors than usual.
If on the other hand, errors of size 1e-6 (say) or more for the direct
methods **DOES** indicate a serious bug.
----------------------------------------------------------------------

A5.  It won't compile on my BrandX machine running X-OS.  What do I do?

Firstly make sure that you run "configure" or "configure --with-..." to get
the options you want.  Then you use "make".  The simplest installation
procedure for Unix systems to create the entire library is

	configure --with-all
	make all