gstat-2.5.1.tar.gz_gstat

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 numb