meschach tutorial.htm

来自「C语言版本的矩阵库」· HTM 代码 · 共 870 行 · 第 1/3 页

    v_linlist(temp,v1,1.0,v2,2.0,v3,2.0,v4,1.0,VNULL)
    /* adjust x */
    v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */

    return t+h;                   /* return the new time */
}
</PRE>
<H2><A name=least_sq>A least squares problem</A></H2>Here we need to use 
matrices and matrix factorisations (in particular, a QR factorisation) in order 
to find the best linear least squares solution to some data. Thus in order to 
solve the (approximate) equations 
<P><I>Ax approx= b,</I>
<P>for <I>x</I> where <I>A</I> is an <I>m x n</I> matrix (<I>m&gt;n</I>) we 
really need to solve the optimisation problem 
<P>min_<I>x ||Ax-b||_2^2.</I> If we write <I>A=QR</I> where <I>Q</I> is an 
orthogonal <I>m x m</I> matrix and <I>R</I> is an upper triangular <I>m x n</I> 
matrix then 
<P><PRE>    ||Ax-b||_2=||Rx-Q^T b||_2 = ||[R_1]b  [Q_1^T]b||
                                ||[ O ] - [Q_2^T] ||_2
</PRE>
<P>where <I>R_1</I> is an <I>n x n</I> upper triangular matrix. If <I>A</I> has 
full rank then <I>R_1</I> will be an invertible matrix, and the best least 
squares solution of <I>Ax approx= b</I> is <I>x=inverse(R_1) Q_1^T b</I>. These 
calculations can be be done quite easily as there is a <TT>QRfactor()</TT> 
function available with the system. <TT>QRfactor()</TT> is declared to have the 
prototype <PRE>MAT     *QRfactor(MAT *A, VEC *diag);
</PRE>The matrix <TT>A</TT> is overwritten with the factorisation of <TT>A</TT> 
``in compact form''; that is, while the upper triangular part of <TT>A</TT> is 
indeed the <I>R</I> matrix described above, the <I>Q</I> matrix is stored as a 
collection of Householder vectors in the strictly lower triangular part of 
<TT>A</TT> and in the <TT>diag</TT> vector. The <TT>QRsolve()</TT> function 
knows and uses this compact form and solves <I>QRx approx= b</I> with the call 
<TT>QRsolve(A,diag,b,x)</TT>, which also returns <TT>x</TT>. Here is the code to 
obtain the matrix <I>A</I>, perform the QR factorisation, obtain the data vector 
<I>b</I>, solve for <I>x</I>, and determine what the norm of the errors 
(<I>||Ax-b||_2</I>) is. <PRE>#include "matrix2.h"

main()
{
    MAT *A, *QR;
    VEC *b, *x, *diag;

    /* read in A matrix */
    printf("Input A matrix:\n");

    A = m_input(MNULL); /* A has whatever size is input */
</PRE><PRE>    if ( A-&gt;m &lt; A-&gt;n )
    {
        printf("Need m &gt;= n to obtain least squares fit\n");
        exit(0);
    }
    printf("# A =\n");       m_output(A);
    diag = v_get(A-&gt;m);
    /* QR is to be the QR factorisation of A */
    QR = m_copy(A,MNULL);
    QRfactor(QR,diag);   
    /* read in b vector */
    printf("Input b vector:\n");
    b = v_get(A-&gt;m);
    b = v_input(b);
    printf("# b =\n");       v_output(b);

    /* solve for x */
    x = QRsolve(QR,diag,b,VNULL);
    printf("Vector of best fit parameters is\n");
    v_output(x);
    /* ... and work out norm of errors... */
    printf("||A*x-b|| = 
           v_norm2(v_sub(mv_mlt(A,x,VNULL),b,VNULL)));
}
</PRE>Note that as well as the usual memory allocation functions like 
<TT>m_get()</TT>, the I/O functions like <TT>m_input()</TT> and 
<TT>m_output()</TT>, and the factorise--and--solve functions <TT>QRfactor()</TT> 
and <TT>QRsolve()</TT>, there are also functions for matrix--vector 
multiplication: <PRE>mv_mlt(MAT *A, VEC *x, VEC *out)</PRE>. and also vector--matrix 
multiplication (with the vector on the left): <PRE>vm_mlt(MAT *A, VEC *x, VEC *out)</PRE>, with <I>out=x^T A</I>. There are 
also functions to perform matrix arithmetic --- matrix addition 
<TT>m_add()</TT>, matrix--scalar multiplication <TT>sm_mlt()</TT>, 
matrix--matrix multiplication <TT>m_mlt()</TT>. Several different sorts of 
matrix factorisation are supported: LU factorisation (also known as Gaussian 
elimination) with partial pivoting, by <TT>LUfactor()</TT> and 
<TT>LUsolve()</TT>. Other factorisation methods include Cholesky factorisation 
<TT>CHfactor()</TT> and <TT>CHsolve()</TT>, and QR factorisation with column 
pivoting <TT>QRCPfactor()</TT>. Pivoting involve <I>permutations</I> which have 
their own <TT>PERM</TT> data structure. Permutations can be created by 
<TT>px_get()</TT>, read and written by <TT>px_input()</TT> and 
<TT>px_output()</TT>, multiplied by <TT>px_mlt()</TT>, inverted by 
<TT>px_inv()</TT> and applied to vectors by <TT>px_vec()</TT>. The above program 
can be put into a file <TT>leastsq.c</TT> and compiled under Unix<I>^{TM}</I> 
using <PRE>cc -o leastsq leastsq.c meschach.a -lm
</PRE>A sample session using <TT>leastsq</TT> follows: <PRE>
Input A matrix:
Matrix: rows cols:5 3
row 0:
entry (0,0): 3
entry (0,1): -1
entry (0,2): 2
Continue: 
row 1:
entry (1,0): 2
entry (1,1): -1
entry (1,2): 1
Continue: n
row 1:
entry (1,0): old              2 new: 2
entry (1,1): old             -1 new: -1
entry (1,2): old              1 new: 1.2
Continue: 
row 2:
entry (2,0): old              0 new: 2.5
  ....
  ....             (Data entry)
  ....
# A =
Matrix: 5 by 3
row 0:              3             -1              2 
row 1:              2             -1            1.2 
row 2:            2.5              1           -1.5 
row 3:              3              1              1 
row 4:             -1              1           -2.2 
Input b vector:
entry 0: old              0 new: 5
entry 1: old              0 new: 3
entry 2: old              0 new: 2
entry 3: old              0 new: 4
entry 4: old              0 new: 6
# b =
Vector: dim: 5
         5          3          2          4          6 
Vector of best fit parameters is
Vector: dim: 3
    1.47241555   -0.402817858    -1.14411815 
||A*x-b|| = 6.78938
</PRE>The <I>Q</I> matrix can be obtained explicitly by the routine 
<TT>makeQ()</TT>. The <I>Q</I> matrix can then be used to obtain an orthogonal 
basis for the range of <I>A</I>. An orthogonal basis for the null space of 
<I>A</I> can be obtained by finding the QR-factorisation of <I>A^T</I>. 
<H2><A name=sparse_eg>A sparse matrix example</A></H2>To illustrate the sparse 
matrix routines, consider the problem of solving Poisson's equation on a square 
using finite differences, and incomplete Cholesky factorisation. The actual 
equations to solve are 
<P><I>u_{i,j+1}+u_{i,j-1}+u_{i+1,j}+u_{i-1,j}-4u_{ij}= h^2 f(x_i,y_j)</I>,
<P>for <I>i,j=1,...,N</I> where <I>u_{0,j}=u_{i,0}=u_{N+1,j}=u_{i,N+1}=0</I> for 
<I>i,j=1,...,N</I> and <I>h</I> is the common distance between grid points. The 
first task is to set up the matrix describing this system of linear equations. 
The next is to set up the right-hand side. The third is to form the incomplete 
Cholesky factorisation of this matrix, and finally to use the sparse matrix 
conjugate gradient routine with the incomplete Cholesky factorisation as 
preconditioner. Setting up the matrix and right-hand side can be done by the 
following code: <PRE>#define N 100
#define index(i,j) (N*((i)-1)+(j)-1)
  ......
A = sp_get(N*N,N*N,5);
b = v_get(N*N);
h = 1.0/(N+1);      /* for a unit square */
  ......
</PRE><PRE>for ( i = 1; i &lt;= N; i++ )
    for ( j = 1; j &lt;= N; j++ )
    {
        if ( i &lt; N )
            sp_set_val(A,index(i,j),index(i+1,j),-1.0);
        if ( i &gt; 1 )
            sp_set_val(A,index(i,j),index(i-1,j),-1.0);
        if ( j &lt; N )
            sp_set_val(A,index(i,j),index(i,j+1),-1.0);
        if ( j &gt; 1 )
            sp_set_val(A,index(i,j),index(i,j-1),-1.0);
        sp_set_val(A,index(i,j),index(i,j),4.0);
        b-&gt;ve[index(i,j)] = -h*h*f(h*i,h*j);
    }
</PRE>Once the matrix and right-hand side are set up, the next task is to 
compute the sparse incomplete Cholesky factorisation of <TT>A</TT>. This must be 
done in a different matrix, so <TT>A</TT> must be copied. <PRE>LLT = sp_copy(A);
spICHfactor(LLT);
</PRE>Now when that is done, the remainder is easy: <PRE>out = v_get(A-&gt;m);
  ......
iter_spcg(A,LLT,b,1e-6,out,1000,&amp;num_steps);
printf("Number of iterations = 
  ......
</PRE>and the output can be used in whatever way desired. For graphical output 
of the results, the solution vector can be copied into a square matrix, which is 
then saved in MATLAB(TM) format using <TT>m_save()</TT>, and graphical output 
can be produced by MATLAB(TM). 
<H2><A name=how_to>How do I ....?</A></H2>For the convenience of the user, here 
a number of common tasks that people need to perform frequently, and how to 
perform the computations using Meschach. 
<H3>....solve a system of linear equations</H3>If you wish to solve <I>Ax=b</I> 
for <I>x</I> given <I>A</I> and <I>b</I> (without destroying <I>A</I>), then the 
following code will do this: <PRE>VEC     *x, *b;
MAT	*A, *LU;
PERM	*pivot;
  ......
LU = m_get(A-&gt;m,A-&gt;n);
LU = m_copy(A,LU);
pivot = px_get(A-&gt;m);
LUfactor(LU,pivot);
/* set values of b here */
x = LUsolve(LU,pivot,b,VNULL);
</PRE>
<H3>....solve a least-squares problem</H3>To minimise 
<I>||Ax-b||_2^2=\sum_i((Ax)_i-b_i)^2</I>, the most reliable method is based on 
the QR-factorisation. The following code performs this calculation assuming that 
<I>A</I> is <I>m</I> x <I>n</I> with <I>m &gt;= n</I>: <PRE>MAT	*A, *QR;
VEC	*diag, *b, *x;
  ......
QR = m_get(A-&gt;m,A-&gt;n);
QR = m_copy(A,QR);
diag = v_get(A-&gt;n);
QRfactor(QR,diag);
/* set values of b here */
x = QRsolve(QR,diag,b,x);
</PRE>
<H3>.... find all the eigenvalues (and eigenvectors) of a general matrix</H3>The 
best method is based on the <I>Schur decomposition</I>. For symmetric matrices, 
the eigenvalues and eigenvectors can be computed by a single call to 
<TT>symmeig()</TT>. For non-symmetric matrices, the situation is more complex 
and the problem of finding eigenvalues and eigenvectors can become quite 
ill-conditioned. Provided the problem is not too ill-conditioned, the following 
code should give accurate results: <PRE>/* A is the matrix whose e-vals and e-vecs are sought */
MAT	*A, *T, *Q, *X_re, *X_im;
VEC	*evals_re, *evals_im;
  ......
Q = m_get(A-&gt;m,A-&gt;n);
T = m_copy(A,MNULL);
/* compute Schur form: A = Q.T.Q^T */
schur(T,Q);
/* extract eigenvalues */
evals_re = v_get(A-&gt;m);
evals_im = v_get(A-&gt;m);
schur_evals(T,evals_re,evals_im);
/* Q not needed for eiegenvalues */
X_re = m_get(A-&gt;m,A-&gt;n);
X_im = m_get(A-&gt;m,A-&gt;n);
schur_vecs(T,Q,X_re,X_im);
/* k'th eigenvector is k'th column of (X_re + i*X_im) */
</PRE>
<H3>.... solve a large, sparse, positive definite system of equations</H3>An 
example of a large, sparse, positive definite matrix is the matrix obtained from 
a finite-difference approximation of the Laplacian operator. If an explicit 
representation of such a matrix is available, then the following code is 
suggested as a reasonable way of computing solutions: <PRE>/* A.x == b is the system to be solved */
sp_mat	*A, *LLT;
VEC	*x, *b;
int     num_steps;
  ......
/* set up A and b */
  ......
x = m_get(A-&gt;m);
LLT = sp_copy(A);
/* preconditioning using incomplete Cholesky */
spICHfactor(LLT);
/* now use pre-conditioned conjugate gradients */
x = iter_spcg(A,LLT,b,1e-7,x,1000,&amp;num_steps);
/* solution computed with relative residual &lt; 10^{-7} */
</PRE>If explicitly storing such a matrix takes up too much memory, then if you 
can write a routine to perform the calculation of <I>Ax</I> for any given 
<I>x</I>, the following code may be more suitable (if slower): <PRE>VEC	*mult_routine(user_def,x,out)
void    *user_def;
VEC	*x, *out;
{
  /* compute out = A*x */
  ......
  return out;
}
</PRE><PRE>main()
{
    ITER *ip;
    VEC  *x, *b;
      ......
    b = v_get(BIG_DIM);  /* right-hand side */
    x = v_get(BIG_DIM);  /* solution */

    /* set up b */
      ......
    ip = iter_get(b-&gt;dim, x-&gt;dim);
    ip-&gt;rhs = v_copy(b,ip-&gt;rhs);
    ip-&gt;info = NULL;     /* if you don't want information
                                about solution process */
    v_zero(ip-&gt;x);       /* initial guess is zero */
    iter_Ax(ip,mult_routine,user_def);
    iter_cg(ip);
    printf("# Solution is:\n");   v_output(ip-&gt;x);
      ......
    ITER_FREE(ip);       /* destroy ip */
}
</PRE>The <TT>user_def</TT> argument is for a pointer to a user-defined 
structure (possibly NULL, if you don't need this) so that you can write a common 
function for handling a large number of different circumstances. 
</LI></BODY></HTML>