meschach tutorial.htm
来自「C语言版本的矩阵库」· HTM 代码 · 共 870 行 · 第 1/3 页
HTM
870 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0070)http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html -->
<HTML><HEAD><TITLE>Meschach: Tutorial</TITLE>
<META http-equiv=Content-Type content="text/html; charset=utf-8">
<META content="MSHTML 6.00.2800.1106" name=GENERATOR></HEAD>
<BODY>
<H1>Tutorial</H1>Sections:
<UL>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#basic_ops">Data
structures and some basic operations</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#mem_manage">How
to manage memory</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#vec_ops">Simple
vector operations: An RK4 routine</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#list_args">Using
routines for lists of arguments</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#least_sq">A
least squares problem</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#sparse_eg">A
sparse matrix example</A>
<LI><A
href="http://www.math.uiowa.edu/~dstewart/meschach/html_manual/tutorial.html#how_to">How
do I ....?</A> </LI></UL>In this chapter, 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 chapter are implemented in
<TT>tutorial.c</TT>. 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.
<H2><A name=basic_ops>The data structures and some basic operations</A></H2>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 <I>pointers</I>
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 <I>A</I> of size <I>3 x 4</I>, a
vector <I>x</I> of dimension 10, and a permutation <I>p</I> of size 10, use the
following code: <PRE>#include "matrix.h"
..............
main()
{
MAT *A;
VEC *x;
PERM *p;
..........
A = m_get(3,4);
x = v_get(10);
p = px_get(10);
..........
}
</PRE>This initialises these data structures to have the given size. The matrix
<I>A</I> and the vector <I>x</I> are initially all zero, while <I>p</I> is
initially the identity permutation. They can be disposed of by calling
<TT>M_FREE(A)</TT>, <TT>V_FREE(x)</TT> and <TT>PX_FREE(p)</TT> 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>(i,j)</I> component of <I>A</I> is
accessed by <TT>A->me[i][j]</TT>, <I>x_i</I> by <TT>x->ve[i]</TT> and
<I>p_i</I> by <TT>p->pe[i]</TT>. Their sizes are also directly accessible:
<TT>A->m</TT> and <TT>A->n</TT> are the number of rows and columns of
<I>A</I> respectively, <TT>x->dim</TT> is the dimension of <I>x</I>, and size
of <I>p</I> is <TT>p->size</TT>. Note that the indexes are <I>zero
relative</I> just as they are in ordinary C, so that the index <TT>i</TT> in
<TT>x->ve[i]</TT> can range from 0 to <TT>x->dim</TT><I>-1</I>. Thus the
total number of entries of a vector is exactly <TT>x->dim</TT>. 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 <I>x</I> to <I>y</I> it is sufficient to write <TT>y =
v_copy(x,VNULL)</TT>. The <TT>VNULL</TT> is the NULL vector, and usually tells
the routine called to create a vector for output. Thus, the <TT>v_copy</TT>
function will create a vector which has the same size as <I>x</I> and all the
components are equal to those of <I>x</I>. If <TT>y</TT> has already been
created then you can write <TT>y = v_copy(x,y)</TT>; in general, writing
``<TT>v_copy(x,y);</TT>'' is not enough! If <TT>y</TT> is NULL, then it is
created (to have the correct size, i.e. the same size as <TT>x</TT>), and if it
is the wrong size, then it is resized to have the correct size (i.e. same size
as <TT>x</TT>). 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:
<TT>v_add(x,y,out)</TT>, <TT>v_sub(x,y,out)</TT>, <TT>sv_mlt(s,x,out)</TT> where
<TT>x</TT> and <TT>y</TT> are input vectors (with data type <TT>VEC *</TT>),
<TT>out</TT> is the output vector (same data type) and <TT>s</TT> is a double
precision number (data type <TT>double</TT>). There is also a special
combination routine, which computes <I>out=v_1+s v_2</I> in a single routine:
<TT>v_mltadd(v1,v2,s,out)</TT>. 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:
<TT>in_prod(x,y)</TT> returns the inner product of <I>x</I> and <I>y</I>. 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: <TT>m_add(A,B,C)</TT> which returns
<I>C=A+B</I>, and <TT>sm_mlt(s,A,C)</TT> which returns <I>C=sA</I>. The data
types of <TT>A</TT>, <TT>B</TT> and <TT>C</TT> are all <TT>MAT *</TT>, while
that of <TT>s</TT> is type <TT>double</TT> as before. The matrix NULL is called
<TT>MNULL</TT>. Multiplying matrices and vectors can be done by a single
function call: <TT>mv_mlt(A,x,out)</TT> returns <I>out=A x</I> while
<TT>vm_mlt(A,x,out)</TT> returns <I>out=A^Tx</I>, or equivalently,
<I>out^T=x^TA</I>. 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 <TT>px_vec(p,x,p_x)</TT>, rows and columns of matrices can be permuted
by using <TT>px_rows(p,A,p_A)</TT>, <TT>px_cols(p,A,A_p)</TT>, and permutations
can be multiplied using <TT>px_mlt(p1,p2,p1_p2)</TT> and inverted using
<TT>px_inv(p,p_inv)</TT>. The NULL permutation is called <TT>PXNULL</TT>. There
are also utility routines to initialise or re-initialise these data structures:
<TT>v_zero(x)</TT>, <TT>m_zero(A)</TT>, <TT>m_ident(A)</TT> (which sets
<I>A=I</I> of the correct size), <TT>v_rand(x)</TT>, <TT>m_rand(A)</TT> which
sets the entries of <TT>x</TT> and <TT>A</TT> respectively to be randomly and
uniformly selected between zero and one, and <TT>px_ident(p)</TT> which sets
<TT>p</TT> to be an identity permutation. Input and output are accomplished by
library routines <TT>v_input(x)</TT>, <TT>m_input(A)</TT>, and
<TT>px_input(p)</TT>. If a null object is passed to any of these input routines,
all data will be obtained from the input file, which is <TT>stdin</TT>. 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:
<TT>v_output(x)</TT>, <TT>m_output(A)</TT> and <TT>px_output(p)</TT>. This
output is both human and machine readable! If you wish to send the data to a
file other than the standard output device <TT>stdout</TT>, or receive input
from a file or device other than the standard input device <TT>stdin</TT>, take
the appropriate routine above, use the ``<TT>foutpout</TT>'' suffix instead of
just ``<TT>output</TT>'', and add a file pointer as the first argument. For
example, to send a matrix <I>A</I> to a file called ``fred'', use the following:
<PRE>#include "matrix.h"
.............
main()
{
FILE *fp;
MAT *A;
.............
fp = fopen("fred","w");
m_foutput(fp,A);
.............
}
</PRE>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: <PRE># The initial vector must not be zero
# x =
Vector: dim: 3
-7 0 3
</PRE>For general input/output which conforms to this format, allowing comments
in the input files, use the <TT>input()</TT> and <TT>finput()</TT> macros. These
are used to print out a prompt message if <TT>stdin</TT> 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: <PRE>input("Input number of steps: ","
fp = stdin;
finput(fp,"Input number of steps: ","
fp = fopen("fred","r");
finput(fp,"Input number of steps: ","
</PRE>The <TT>"%d"</TT> strings are the format strings as used by
<TT>scanf()</TT> and <TT>fscanf()</TT>; the last argument is the pointer to the
variable (unless the variable is a string) just as for <TT>scanf()</TT> and
<TT>fscanf()</TT>. The first two macro calls read input from <TT>stdin</TT>, the
last from the file <TT>fred</TT>. If, in the first two calls, <TT>stdin</TT> is
a keyboard (a ``tty'' in Unix jargon) then the prompt string <TT>"Input number
of steps: "</TT> is printed out on the terminal. The second part of the library
contains routines for various factorisation methods. To use it put <PRE>#include "matrix2.h"
</PRE>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 <I>Q</I> 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 <I>LDL^T</I> factorisations. For using complex numbers,
vectors and matrices include <PRE>#include "zmatrix.h"
</PRE>for using the basic routines, and <PRE>#include "zmatrix2.h"
</PRE>for the complex matrix factorisation routines. The <TT>zmatrix2.h</TT>
file includes <TT>matrix.h</TT> and <TT>zmatrix.h</TT> so you don't need these
files included together. For using the sparse matrix routines in the library you
need to put <PRE>#include "sparse.h"
</PRE>or, if you use any sparse factorisation routines <PRE>#include "sparse2.h"
</PRE>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 <PRE>#include "iter.h"
</PRE>This includes the <TT>sparse.h</TT> and <TT>matrix.h</TT> 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
(<TT>sin()</TT>, <TT>cos()</TT>, <TT>tan()</TT>, <TT>exp()</TT>, <TT>log()</TT>,
<TT>sqrt()</TT>, <TT>acos()</TT> etc.) don't forget to include the standard
mathematics header: <PRE>#include <MATH.H>
</PRE>This file is <I>not</I> automatically included by any of the Meschach
header files.
<H2><A name=mem_manage>How to manage memory</A></H2>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.
<P><NL>
<LI>Memory allocation, deallocation and resizing takes a significant amount of
time compared with (say) vector operations, so it should not be done too
frequently.
<LI>Allocating memory but not deallocating it means that it can't 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. </NL>
<P>There are three main strategies that are recommended for deciding how to
allocate, deallocate and resize objects. These are ``<I>no deallocation</I>''
which is really only useful for demonstration programs, ``<I>allocate and
deallocate</I>'' which minimises overall memory requirements at the expense of
speed, and ``<I>resize on demand</I>'' which is useful for routines that are
called repeatedly. A new technique for static workspace arrays is to
``<I>register workspace variables</I>''.
<H3>No deallocation</H3>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 <PRE>QR = m_copy(A,MNULL); /* allocate memory for QR */
</PRE>to allocate the memory, but without the call <TT>M_FREE(QR);</TT> in it.
This can be acceptable if <TT>QR = m_copy(A,MNULL)</TT> is only executed once,
and so the allocated memory never needs to be explicitly deallocated. This would
<I>not</I> be acceptable if <TT>QR = m_copy(A,MNULL)</TT> occurred inside a
<TT>for</TT> 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 <TT>QR</TT>. The next subsection shows how to avoid this.
<H3>Allocate and deallocate</H3>This is the most straightforward way of ensuring
that memory is not lost. With the example of allocating <TT>QR</TT> it would
work like this: <PRE>for ( ... ; ... ; ... )
{
QR = m_copy(A,MNULL); /* allocate memory for QR */
/* could have been allocated by m_get() */
/* use QR */
......
......
/* deallocate QR so memory can be reused */
M_FREE(QR);
}
</PRE>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).
<H3>Resize on demand</H3>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 <TT>v1</TT>,
<TT>v2</TT>, etc. in the Runge--Kutta routine <TT>rk4()</TT> are allocated
according to this strategy: <PRE>rk4(...,x,...)
{
static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL,
*v4=VNULL, *temp=VNULL;
.......
v1 = v_resize(v1,x->dim);
v2 = v_resize(v2,x->dim);
v3 = v_resize(v3,x->dim);
v4 = v_resize(v4,x->dim);
temp = v_resize(temp,x->dim);
.......
}
</PRE>The intention is that the <TT>rk4()</TT> routine is called repeatedly with
the same size <TT>x</TT> vector. It then doesn't make as much sense to allocate
<TT>v1</TT>, <TT>v2</TT> etc. whenever the function is called. Instead,
<TT>v_resize()</TT> only performs memory allocation if the memory already
allocated to <TT>v1</TT>, <TT>v2</TT> etc. is smaller than <TT>x->dim</TT>.
The vectors <TT>v1</TT>, <TT>v2</TT> etc. are declared to be <TT>static</TT> to
ensure that their values are not lost between function calls. Variables that are
declared <TT>static</TT> are set to <TT>NULL</TT> or zero by default. So the
declaration of <TT>v1</TT>, <TT>v2</TT>, etc., could be <PRE>static VEC *v1, *v2, *v3, *v4, *temp;
</PRE>This strategy of resizing static workspace variables is not so useful if
the object being allocated is extremely large. The previous ``allocate and
deallocate'' strategy is much more efficient for memory in those circumstances.
However, the following section shows how to get the best of both worlds.
<H3>Registration of workspace</H3>From version 1.2 onwards, workspace variables
can be <I>registered</I> so that the memory they reference can be freed up on
demand. To do this, the function containing the static workspace variables has
to include calls to <TT>MEM_STAT_REG(var,type)</TT> where <TT>var</TT> is a
pointer to a Meschach data type (such as <TT>VEC</TT> or <TT>MAT</TT>). This
call should be placed <I>after</I> the call to the appropriate resize function.
The <TT>type</TT> parameter should be a <TT>TYPE_...</TT> macro where the
``<TT>...</TT>'' is the name of a Meschach type such as <TT>VEC</TT> or
<TT>MAT</TT>. For example, <PRE>rk4(...,x,...)
{
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?