faq_meschach.txt

C语言版本的矩阵库

If you want to use gcc (or if that is not the default compiler) then use

	configgnu --with-all
	make all

The "configure" script is not perfect and sometimes misses some things.
For example, when we used configure on Linux, it missed the fact that the
type u_int was already declared.  This meant that some editing of the
generated machine.h file was needed.

Users of MS-DOS machines and C compilers, and users of VMS and other
non-Unix operating systems will need to manually edit your machine.h and/or
makefile for your system.  Some guidance in this can come from the examples
in the MACHINES/xxxx directories which contain machine.h, makefile (and
occasionally machine.c) files for some machines.

If configure fails, then the files to modify for your particular machine
are "machine.h" and "makefile".  (The file "machine.c" can also be
modified for speed.)  There are a number of macros in machine.h which can
be set to control the way the library is compiled.  You should consult a
local expert (and/or your local manuals!) to determine how things should be
set.  If this is not successful, send us a message.  We are unlikely to
know the magic words for your particular machine/compiler/operating system,
but we may know someone who does.
----------------------------------------------------------------------

A6.  Is there a C++ version of Meschach?

There is a C++ version (called Meschach++) which is under development.  
A beta version can be found in the directory MESCHPP.
There is a manual for this version.  Note that this manual will not 
be a substitute for the main manual described above -- the Meschach++ 
manual will refer to "Meschach: Matrix Computations in C" for many 
details on the behaviour of the various routines.  The C++ version is, 
basically, a "wrapper" for Meschach.  

Meschach++ is not being developed by David Stewart and Zbigniew Leyk, but
by Stephen Roberts, Brian Austin and Alex Austin also while at the
Australian National University, School of Mathematical Sciences.
You can contact Stephen Roberts at [e-mail].

Some points should be made about the use of C++.  The use of operator
overloading is one of the most attractive features of C++, especially for
tasks such as matrix computations.  However, it also has some drawbacks.
One drawback is the problem of temporary variables.  In a matrix expression
of the form:

	S = A + B + C + D + E;

where all the objects are matrices (of the same size) if operator
overloading is implemented in the naive way, then this will result in
temporary matrices for (A+B), ((A+B)+C), (((A+B)+C)+D) and
((((A+B)+C)+D)+E).  In fact, the sum could be accumulated directly into S
(provided there is no aliasing) and there do not need to be any
temporaries.  Instead, a naive implementation would spend most of its time
allocating and deallocating temporary variables.  This is just one aspect
of a problem with optimising evaluation of matrix expressions.

There are ways of fixing this problem using delayed evaluation, which can
be implemented in C++.  .  Another is to use 3-term functions, similar in
spirit to Meschach's use of a final output parameter.  This loses the
elegance of operator overloading, but gains greatly in performance.
----------------------------------------------------------------------

A7.  Meschach doesn't have an xxxx routine.  Why not?  What do I do?

There is usually one answer to "Why not?"  We haven't had time.
But some operations should be avoided if possible.  It was a long time
before we implemented matrix inversion as a single function.  Most of the
time it is not needed: often only the solution of systems of equations is
needed.  There is no single "eigenvalue/vector" routine.  This is because
this is done by computing the Schur decomposition of a matrix, and then
extracting information about the eigenvalues and eigenvectors.  The Schur
decomposition not only provides more information, but in the case of
eigenvectors, is much more numerically stable.  After all, not all matrices
can be transformed into diagonal form -- they have fewer eigenvectors than
columns. 

The other point that should be made is that Meschach is a toolkit.
Often it is up to you to make the part you need.  You may need to
re-consider what you are doing -- is it numerically stable?
Is it an efficient way of doing the task?  It may be better to build on
what Meschach provides, rather than hoping that it will all be there.
----------------------------------------------------------------------

A8.  Can I use Meschach with MATLAB?

The answer is both "Yes" and "No".  The "Yes" part is that Meschach
provides input and output of dense matrices and vectors that is compatible
with MATLAB.  This means that matrices can be transferred between MATLAB
and Meschach programs.  However, Meschach has not been made into a
MEX-compatible library.  So Meschach routines are not **yet** ready to be
used in MEX routines.
----------------------------------------------------------------------

A9.  I don't want ALL of Meschach.  How do I get just a subset?

Yes you can.  If you have the shar files from netlib then the most basic
library can be built out of shar files meschach0.shar and meschach1.shar.
This contains the core routines and the basic matrix, vector computations.
(It is not even enough to get torture to work.)  The next shar file needed
is meschach2.shar which contains the matrix factorisation routines.
With these three shar files the basic library can be built.  The shar file
meschach3.shar contains the sparse and iterative methods, while
meschach4.shar contains the complex vector and matrix operations.

If you want the basic library (with the matrix factorisations), get the
shar files meschach[012].shar, expand them, and then build the library
using

	configure
	make basic

If you want to add just the sparse routines use

	configure --with-sparse
	make sparse

If you want to add just the complex routines use

	configure --with-complex
	make complex
----------------------------------------------------------------------

A10. What does Meschach mean?

Meschach doesn't stand for anything in particular.  And it is NOT an
acronym.  It is, however, someone's name.  If you want to find out, look
two thirds down page ii of the manual :-)
----------------------------------------------------------------------

A11. My code that uses Meschach seems to run very slowly/uses incredible
	amounts of memory.  Why?

It is quite possible to make Meschach run very inefficiently.  What it is
most likely doing is spending all of its time allocating matrices.
(If it is not also deallocating them **explicitly**, then you will use an
amazing amount of memory.)  Look at where you call Meschach.  Is the output
parameter NULL?  If so, then a new data structure is created to hold the
output whenever the routine is called.

For example,

	MAT  *A, *B, *C;

	for ( i = 0; i < 100; i++ )
	{
	    ....
	    C = m_add(A,B,MNULL);	/* C = A+B */
	    ....
	}

will create a new C every time through the loop.  To avoid this use C as
the output parameter:

	MAT  *A, *B, *C;

	C = m_get(A->m,A->n);
	for ( i = 0; i < 100; i++ )
	{
	    ....
	    C = m_add(A,B,C);
	    ....
	}

See the pages under the index entry "memory management" for a comprehensive
discussion of ways of dealing with memory allocation and deallocation.
----------------------------------------------------------------------

A12. I have a huge problem and I keep running out of memory.  What do I do?

Large problems often require specialised techniques.  A large problem is
one where just using dense matrices is intolerably slow.  For most
workstations and most problems this is somewhere around 1000 x 1000,
possibly a bit less, possibly a bit more.

The first strategy to use is to consider using sparse matrices.  Direct
methods can be used with sparse matrices, at least for solving systems of
equations.  The danger here is that fill-in will result in the sparse
matrix becoming a bit too dense for the amount of memory available.
Re-ordering techniques should be used (unfortunately Meschach does not
currently implement these) to reduce the amount of fill-in created by the
factorisation process.

If this is not enough, then iterative methods should be used.  The ones
implemented in Meschach are based on "Krylov subspaces".  These just
require a function to compute the product A.x for any given vector x
(or sometimes the transposed-matrix-vector product A^T.x).  These do not
require the matrix involved to be sparse; it can be dense, but as long
as it has the right structure so that A.x can be computed quickly,
these methods are appropriate.  Preconditioning is desirable for these
iterative methods, and can be achieved in a variety of ways.  One is to
use incomplete factorisations (which simply ignore fill-in).  Other methods
are problem dependent: for partial differential equations, multigrid
preconditioners can work very well.  Classical methods (Gauss-Seidel, SOR,
SSOR, Chebyshev acceleration) can be used to create pre-conditioners which
when combined with a Krylov subspace method can be far better than either
alone.
----------------------------------------------------------------------

A13. How did you come to write Meschach?

The first of us to start working on Meschach was David Stewart.  The first
ideas of how to use pointers and dynamic memory allocation for storing
matrices was developing while completing an Electrical Engineering degree.
The next year was spent as a systems programmer for 5 months and then
starting an Honours year in Mathematics (in USA read "Masters", in Germany
read "Diplom").  By 1986, when he started his PhD he was ready to start
writing code for it.  The library was never an end in itself, but a tool
for writing the real programs for the thesis (in numerical analysis and
differential equations).

In 1989 the library was the basis for programs using linear programming and
ODE solvers; however, it was still very limited and had no
eigenvalue/vector routines, no complex routines, and was not particularly
portable.  For a while he was waiting for a job on a project at an
operations research company on large scale optimisation -- the project
eventually fell through, but in the mean time some basic sparse routines
were written in anticipation.  Later that year David Stewart was hired as a
lecturer in Computational Mathematics at the University of Queensland,
Australia, and taught Numerical Linear Algebra.  As a result, the
development of a number of routines took place in parallel with them being
taught in the classroom.

Meanwhile, Zbigniew Leyk had also completed his PhD and was working in a
Post-doc at Cornell.  While there he developed his own ideas about a
numerical library in C as part of his work on multigrid methods.  This had
a different flavour to Meschach in some respects, but both approaches
emphasised the use of data structures.

By 1991 when both David Stewart and Zbigniew Leyk arrived at the Australian
National University, the ideas had matured enough; Meschach had been
developed in terms of scope and portability to become a really general
tool.  Meschach was placed in netlib in 1991.  A manual had been written
over 1990/1991.  There was mutual interest in working with Meschach from
both of us, and during 1993 we worked intensively on a major new revision,
which has become version 1.2.  This had complex data structures and
operations, a revamped way of dealing with iterative methods, better
methods for handling workspace vectors and matrices and a new installation
procedure based on GNU's autoconf.
----------------------------------------------------------------------

A14. Complex part of Meschach

>> I wrote for myself a file called complex.h
>> and this is the one which is included 
>> automatically in your files. How can this work?
>> What type of complex.h is required? 

complex.h is not for a C++ header file.  Rather, the idea is that if a user
have a better local complex.h (C) header file then that should be used.
However, there may be problems with names being different (e.g.  cadd() vs.
zadd() -- Meschach uses the latter).  It might be necessary to edit
machine.h to make it work better for a particular setup.
--------------------------------------------------------------------------

A15. Can I do an index search in order to know what routines exist?

In meschach0.shar (and in the tar and zip distributions) is a file called
fnindex.txt which contains a list of routines and briefly what they do.
Use your favourite string search method to find the right routine(s).
------------------------------------------------------------------------

End of Meschach FAQ