tutorial.txt

C语言版本的矩阵库

			   MESCHACH VERSION 1.2A
			   ---------------------


				 TUTORIAL
				 ========


   In this manual the basic data structures are introduced, and some of the
more basic operations are illustrated.  Then some examples of how to use
the data structures and procedures to solve some simple problems are given.
The first example program is a simple 4th order Runge-Kutta solver for
ordinary differential equations.  The second is a general least squares
equation solver for over-determined equations.  The third example
illustrates how to solve a problem involving sparse matrices.  These
examples illustrate the use of matrices, matrix factorisations and solving
systems of linear equations.  The examples described in this manual are
implemented in tutorial.c.

   While the description of each aspect of the system is brief and far from
comprehensive, the aim is to show the different aspects of how to set up
programs and routines and how these work in practice, which includes I/O
and error-handling issues.



1.  THE DATA STRUCTURES AND SOME BASIC OPERATIONS

   The three main data structures are those describing vectors, matrices
and permutations.  These have been used to create data structures for
simplex tableaus for linear programming, and used with data structures for
sparse matrices etc.  To use the system reliably, you should always use
pointers to these data structures and use library routines to do all the
necessary initialisation.

   In fact, for the operations that involve memory management (creation,
destruction and resizing), it is essential that you use the routines
provided.

   For example, to create a matrix A of size 34 , a vector x of dimension
10, and a permutation p of size 10, use the following code:


  #include "matrix.h"
  ..............
  main()
  {
     MAT   *A;
     VEC   *x;
     PERM  *p;
     ..........
     A = m_get(3,4);
     x = v_get(10);
     p = px_get(10);
     ..........
  }


   This initialises these data structures to have the given size.  The
matrix A and the vector x are initially all zero, while p is initially the
identity permutation.

   They can be disposed of by calling M_FREE(A), V_FREE(x) and PX_FREE(p)
respectively if you need to re-use the memory for something else.  The
elements of each data structure can be accessed directly using the members
(or fields) of the corresponding structures.  For example the (i,j)
component of A is accessed by A->me[i][j], x_i by x->ve[i] and p_i by
p->pe[i].

   Their sizes are also directly accessible: A->m and A->n are the number
of rows and columns of A respectively, x->dim is the dimension of x , and
size of p is p->size.

   Note that the indexes are zero relative just as they are in ordinary C,
so that the index i in x->ve[i] can range from 0 to x->dim -1 .  Thus the
total number of entries of a vector is exactly x->dim.

   While this alone is sufficient to allow a programmer to do any desired
operation with vectors and matrices it is neither convenient for the
programmer, nor efficient use of the CPU.  A whole library has been
implemented to reduce the burden on the programmer in implementing
algorithms with vectors and matrices.  For instance, to copy a vector from
x to y it is sufficient to write y = v_copy(x,VNULL).  The VNULL is the
NULL vector, and usually tells the routine called to create a vector for
output.

   Thus, the v_copy function will create a vector which has the same size
as x and all the components are equal to those of x.  If y has already
been created then you can write y = v_copy(x,y); in general, writing
``v_copy(x,y);'' is not enough!  If y is NULL, then it is created (to have
the correct size, i.e. the same size as x), and if it is the wrong size,
then it is resized to have the correct size (i.e. same size as x).  Note
that for all the following functions, the output value is returned, even if
you have a non-NULL value as the output argument.  This is the standard
across the entire library.

   Addition, subtraction and scalar multiples of vectors can be computed by
calls to library routines: v_add(x,y,out), v_sub(x,y,out), sv_mlt(s,x,out)
where x and y are input vectors (with data type VEC *), out is the output
vector (same data type) and s is a double precision number (data type
double).  There is also a special combination routine, which computes
out=v_1+s,v_2 in a single routine: v_mltadd(v1,v2,s,out).  This is not only
extremely useful, it is also more efficient than using the scalar-vector
multiply and vector addition routines separately.

   Inner products can be computed directly: in_prod(x,y) returns the inner
product of x and y.  Note that extended precision evaluation is not
guaranteed.  The standard installation options uses double precision
operations throughout the library.

   Equivalent operations can be performed on matrices: m_add(A,B,C) which
returns C=A+B , and sm_mlt(s,A,C) which returns C=sA .  The data types of
A, B and C are all MAT *, while that of s is type double as before.  The
matrix NULL is called MNULL.

   Multiplying matrices and vectors can be done by a single function call:
mv_mlt(A,x,out) returns out=A*x while vm_mlt(A,x,out) returns out=A^T*x , or
equivalently, out^T=x^T*A .  Note that there is no distinction between row
and column vectors unlike certain interactive environments such as MATLAB
or MATCALC.

   Permutations are also an essential part of the package.  Vectors can be
permuted by using px_vec(p,x,p_x), rows and columns of matrices can be
permuted by using px_rows(p,A,p_A), px_cols(p,A,A_p), and permutations can
be multiplied using px_mlt(p1,p2,p1_p2) and inverted using px_inv(p,p_inv).
The NULL permutation is called PXNULL.

   There are also utility routines to initialise or re-initialise these
data structures: v_zero(x), m_zero(A), m_ident(A) (which sets A=I of the
correct size), v_rand(x), m_rand(A) which sets the entries of x and A
respectively to be randomly and uniformly selected between zero and one,
and px_ident(p) which sets p to be an identity permutation.

   Input and output are accomplished by library routines v_input(x),
m_input(A), and px_input(p).  If a null object is passed to any of these
input routines, all data will be obtained from the input file, which is
stdin.  If input is taken from a keyboard then the user will be prompted
for all the data items needed; if input is taken from a file, then the
input will have to be of the same format as that produced by the output
routines, which are: v_output(x), m_output(A) and px_output(p).  This
output is both human and machine readable!

   If you wish to send the data to a file other than the standard output
device stdout, or receive input from a file or device other than the
standard input device stdin, take the appropriate routine above, use the
``foutpout'' suffix instead of just ``output'', and add a file pointer as
the first argument.  For example, to send a matrix A to a file called
``fred'', use the following:


  #include   "matrix.h"
  .............
  main()
  {
     FILE  *fp;
     MAT   *A;
     .............
     fp = fopen("fred","w");
     m_foutput(fp,A);
     .............
  }


   These input routines allow for the presence of comments in the data.  A
comment in the input starts with a ``hash'' character ``#'', and continues
to the end of the line.  For example, the following is valid input for a
3-dimensional vector:

  # The initial vector must not be zero
  # x =
  Vector: dim: 3
  -7      0     3


   For general input/output which conforms to this format, allowing
comments in the input files, use the input() and finput() macros.  These
are used to print out a prompt message if stdin is a terminal (or ``tty''
in Unix jargon), and to skip over any comments if input is from a
non-interactive device.  An example of the usage of these macros is:

  input("Input number of steps: ","%d",&steps);
  fp = stdin;
  finput(fp,"Input number of steps: ","%d",&steps);
  fp = fopen("fred","r");
  finput(fp,"Input number of steps: ","%d",&steps);

The "%d" is one of the format specifiers which are used in fscanf(); the
last argument is the pointer to the variable (unless the variable is a
string) just as for scanf() and fscanf().  The first two macro calls read
input from stdin, the last from the file fred.  If, in the first two calls,
stdin is a keyboard (a ``tty'' in Unix jargon) then the prompt string
  "Input number of steps: " 
is printed out on the terminal.


   The second part of the library contains routines for various
factorisation methods.  To use it put

  #include   "matrix2.h"

at the beginning of your program.  It contains factorisation and solution
routines for LU, Cholesky and QR-factorisation methods, as well as update
routines for Cholesky and QR factorisations.  Supporting these are a number
of Householder transformation and Givens' rotation routines.  Also there is
a routine for generating the Q matrix for a QR-factorisation, if it is
needed explicitly, as it often is.
There are routines for band factorisation and solution for LU and  LDL^T 
factorisations.

For using complex numbers, vectors and matrices include

  #include   "zmatrix.h"

for using the basic routines, and

  #include   "zmatrix2.h"

for the complex matrix factorisation routines.  The zmatrix2.h file
includes matrix.h and zmatrix.h so you don't need these files included
together.

For using the sparse matrix routines in the library you need to put

  #include   "sparse.h"

or, if you use any sparse factorisation routines,

  #include   "sparse2.h"

at the beginning of your file.  The routines contained in the library
include routines for creating, destroying, initialising and updating sparse
matrices, and also routines for sparse matrix-dense vector multiplication,
sparse LU factorisation and sparse Cholesky factorisation.

For using the iterative routines you need to use

  #include   "iter.h"

This includes the sparse.h and matrix.h file.
There are also routines for applying iterative methods such as
pre-conditioned conjugate gradient methods to sparse matrices.

   And if you use the standard maths library (sin(), cos(), tan(), exp(),
log(), sqrt(), acos() etc.)  don't forget to include the standard
mathematics header:

  #include  <math.h>

This file is  not  automatically included by any of the Meschach
header files.



2.  HOW TO MANAGE MEMORY

   Unlike many other numerical libraries, Meschach allows you to allocate,
deallocate and resize the vectors, matrices and permutations that you are
using.  To gain maximum benefit from this it is sometimes necessary to
think a little about where memory is allocated and deallocated.  There are
two reasons for this.

   Memory allocation, deallocation and resizing takes a significant amount
of time compared with (say) vector operations, so it should not be done too
frequently.  Allocating memory but not deallocating it means that it cannot
be used by any other data structure.  Data structures that are no longer
needed should be explicitly deallocated, or kept as static variables for
later use.  Unlike other interpreted systems (such as Lisp) there is no
implicit ``garbage collection'' of no-longer-used memory.

   There are three main strategies that are recommended for deciding how to
allocate, deallocate and resize objects.  These are ``no deallocation''
which is really only useful for demonstration programs, ``allocate and
deallocate'' which minimises overall memory requirements at the expense of
speed, and ``resize on demand'' which is useful for routines that are
called repeatedly.  A new technique for static workspace arrays is to
``register workspace variables''.


2.1  NO DEALLOCATION

   This is the strategy of allocating but never deallocating data
structures.  This is only useful for demonstration programs run with small
to medium size data structures.  For example, there could be a line

  QR = m_copy(A,MNULL);     /* allocate memory for QR */

to allocate the memory, but without the call M_FREE(QR); in it.  This can
be acceptable if QR = m_copy(A,MNULL) is only executed once, and so the
allocated memory never needs to be explicitly deallocated.

   This would not be acceptable if QR = m_copy(A,MNULL) occurred inside a
for loop.  If this were so, then memory would be ``lost'' as far as the
program is concerned until there was insufficient space for allocating the
next matrix for QR.  The next subsection shows how to avoid this.


2.2  ALLOCATE AND DEALLOCATE

   This is the most straightforward way of ensuring that memory is not
lost.  With the example of allocating QR it would work like this:

  for ( ... ; ... ; ... )
  {
    QR = m_copy(A,MNULL);    /* allocate memory for QR */
                             /* could have been allocated by m_get() */
    /* use QR */
      ......
      ......
    /* no longer need QR for this cycle */
    M_FREE(QR);             /* deallocate QR so memory can be reused */
  }

   The allocate and deallocate statements could also have come at the
beginning and end of a function or procedure, so that when the function
returns, all the memory that the function has allocated has been
deallocated.

   This is most suitable for functions or sections of code that are called
repeatedly but involve fairly extensive calculations (at least a
matrix-matrix multiply, or solving a system of equations).


2.3  RESIZE ON DEMAND

   This technique reduces the time involved in memory allocation for code
that is repeatedly called or used, especially where the same size matrix or
vector is needed.  For example, the vectors v1, v2, etc. in the
Runge-Kutta routine rk4() are allocated according to this strategy: