tutorial.txt

C语言版本的矩阵库

  MAT	*A, *LU;
  PERM	*pivot;
  ......
  LU = m_get(A->m,A->n);
  LU = m_copy(A,LU);
  pivot = px_get(A->m);
  LUfactor(LU,pivot);
  /* set values of b here */
  x = LUsolve(LU,pivot,b,VNULL);


7.2  .... SOLVE A LEAST-SQUARES PROBLEM ?

   To minimise ||Ax-b||_2^2 = sum_i ((Ax)_i-b_i)^2, the most reliable
method is based on the QR-factorisation.  The following code performs this
calculation assuming that A is m x n with m > n :

  MAT	*A, *QR;
  VEC	*diag, *b, *x;
  ......
  QR = m_get(A->m,A->n);
  QR = m_copy(A,QR);
  diag = v_get(A->n);
  QRfactor(QR,diag);
  /* set values of b here */
  x = QRsolve(QR,diag,b,x);


7.3  .... FIND ALL THE EIGENVALUES (AND EIGENVECTORS) OF A GENERAL MATRIX ?

   The best method is based on the Schur decomposition.  For symmetric
matrices, the eigenvalues and eigenvectors can be computed by a single call
to symmeig().  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:


  /* A is the matrix whose eigenvalues and eigenvectors are sought */
  MAT	*A, *T, *Q, *X_re, *X_im;
  VEC	*evals_re, *evals_im;
  ......
  Q = m_get(A->m,A->n);
  T = m_copy(A,MNULL);

  /* compute Schur form: A = Q*T*Q^T */
  schur(T,Q);
 
  /* extract eigenvalues */
  evals_re = v_get(A->m);
  evals_im = v_get(A->m);
  schur_evals(T,evals_re,evals_im);

  /* Q not needed for eiegenvalues */
  X_re = m_get(A->m,A->n);
  X_im = m_get(A->m,A->n);
  schur_vecs(T,Q,X_re,X_im);
  /* k'th eigenvector is k'th column of (X_re + i*X_im) */



7.4  .... SOLVE A LARGE, SPARSE, POSITIVE DEFINITE SYSTEM OF EQUATIONS ?

   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:


  /* A*x == b is the system to be solved */
  SPMAT *A, *LLT;
  VEC	*x, *b;
  int   num_steps;
  ......
  /* set up A and b */
  ......
  x = m_get(A->m);
  LLT = sp_copy(A);

  /* preconditioning using the incomplete Cholesky factorisation */
  spICHfactor(LLT);

  /* now use pre-conditioned conjugate gradients */
  x = iter_spcg(A,LLT,b,1e-7,x,1000,&num_steps);
  /* solution computed to give a relative residual of 10^-7 */


   If explicitly storing such a matrix takes up too much memory, then if
you can write a routine to perform the calculation of A*x for any given x ,
the following code may be more suitable (if slower):


  VEC  *mult_routine(user_def,x,out)
  void *user_def;
  VEC  *x, *out;
  {
     /* compute out = A*x */
     ......
     return out;
  }


  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->dim, x->dim);
    ip->b = v_copy(b,ip->b);
    ip->info = NULL;        /* if you don't want information
                                   about solution process */
    v_zero(ip->x);          /* initial guess is zero */
    iter_Ax(ip,mult_routine,user_def);
    iter_cg(ip);
    printf("# Solution is:\n");   v_output(ip->x);
      ......
    ITER_FREE(ip);          /* destroy ip */
  }

   The user_def 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.



8. MORE ADVANCED TOPICS

   Read this if you are interested in using Meschach library as a base for
applications. As an example we show how to implement a new type for 3
dimensional matrices and incorporate this new type into the Meschach
system. Usually this part of Meschach is transparent to a user.  But a more
advanced user can take advantage of these routines. We do not describe
the routines in detail here, but we want to give a rather broad picture of
what can be done.  By the system we mainly mean the system of delivering
information on the number of bytes of allocated memory and routines for
deallocating static variables by mem_stat_... routines.

   First we introduce a concept of a list of types. By a list of types we
mean a set of different types with corresponding routines for creating
these types, destroying and resizing them.  Each type list has a number.
The list 0 is a list of standard Meschach types such as MAT or VEC. Other
lists can be defined by a user or a application (based on Meschach). The
user can attach his/her own list to the system by the routine
mem_attach_list(). Sometimes it is worth checking if a list number is
already used by another application. It can be done by
mem_is_list_attached(ls_num), which returns TRUE if the number ls_num 
is used. And such a list can be removed from the system by
mem_free_list(ls_num) if necessary.

   We describe arguments required by mem_attach_list(). The prototype of
this function is as follow
  
 int mem_attach_list(int ls_num, int ntypes, char *type_names[],
	             int (*free_funcs[])(), MEM_ARRAY sum[]);

where the structure MEM_ARRAY has two members: "bytes" of type long and
"numvar" of type int.  The frst argument is the list number.  Note that you
cannot overwrite another list.  To do this remove first the old list (by
mem_free_list()) or choose another number.  The next argument is the number
of types which are on the list.  This number cannot be changed during
running a program. The third argument is an array containing the names of
types (these are character strings).  The fourth one is an array of
functions deallocating variables of the corresponding type.  And the last
argument is the local array where information about the number of bytes of
allocated/deallocated memory (member bytes) and the number of allocated
variables (member numvar) are gathered. The functions which send
information to this array are mem_bytes_list() and mem_numvar_list().


Example:  The routines described here are in the file tutadv.c.
Firstly we define some macros and a type for 3 dimensional matrices.

#include "matrix.h"
#define M3D_LIST    3      /* list number */
#define TYPE_MAT3D  0      /* the number of a type */
/* type for 3 dimensional matrices */
typedef struct {
	int l,m,n;    /* actual dimensions */
	int max_l, max_m, max_n;    /* maximal dimensions */
	Real ***me;    /* pointer to matrix elements */
	               /* we do not consider segmented memory */
        Real *base, **me2d;  /* me and me2d are additional pointers 
				to base */
} MAT3D;


Now we need two routines: one for allocating memory for 3 dimensional
matrices and the other for deallocating it. It can be useful to have a
routine for resizing 3 dimensional matrices but we do not use it here.
Note the use of mem_bytes_list() and mem_numvar_list() to notify the change
in the number of structures and bytes in use.

/* function for creating a variable of MAT3D type */

MAT3D *m3d_get(l,m,n)
int l,m,n;
{
  MAT3D *mat;
  ....
  /* alocate memory for structure */
  if ((mat = NEW(MAT3D)) == (MAT3D *)NULL)
    error(E_MEM,"m3d_get");
  else if (mem_info_is_on()) {
	/* record how many bytes are allocated to structure */
    mem_bytes_list(TYPE_MAT3D,0,sizeof(MAT3D),M3D_LIST);
	/* record a new allocated variable */
    mem_numvar_list(TYPE_MAT3D,1,M3D_LIST);
  }
  ....
  /* allocate memory for 3D array */
  if ((mat->base = NEW_A(l*m*n,Real)) == (Real *)NULL) 
    error(E_MEM,"m3d_get");
  else if (mem_info_is_on())
    mem_bytes_list(TYPE_MAT3D,0,l*m*n*sizeof(Real),M3D_LIST);
  ....
  return mat;
}

/* deallocate a variable of type MAT3D */

int m3d_free(mat)
MAT3D *mat;
{
  /* do not try to deallocate the NULL pointer */
  if (mat == (MAT3D *)NULL)
    return -1;
  ....
  /* first deallocate base */
  if (mat->base != (Real *)NULL) {
    if (mem_info_is_on())
	/* record how many bytes is deallocated */
      mem_bytes_list(TYPE_MAT3D,mat->max_l*mat->max_m*mat->max_n*sizeof(Real),
		     0,M3D_LIST);
    free((char *)mat->base);
  }
  ....
  /* deallocate  MAT3D structure */
  if (mem_info_is_on()) {
    mem_bytes_list(TYPE_MAT3D,sizeof(MAT3D),0,M3D_LIST);
    mem_numvar_list(TYPE_MAT3D,-1,M3D_LIST);
  }
  free((char *)mat);

  ....
  free((char *)mat);

  return 0;
}


We can now create the arrays necessary for mem_attach_list(). Note that
m3d_sum can be static if it is in the same file as main(), where
mem_attach_list is called. Otherwise it must be global.


char *m3d_names[] = {
  "MAT3D"
};

#define M3D_NUM  (sizeof(m3d_names)/sizeof(*m3d_names))

int (*m3d_free_funcs[M3D_NUM])() = {
  m3d_free
}

static MEM_ARRAY m3d_sum[M3D_NUM];


The last thing is to attach the list to the system.

void main()
{
  MAT3D *M;
  ....

  mem_info_on(TRUE);    /* switch memory info on */
  /* attach the new list */
  mem_attach_list(M3D_LIST,M3D_NUM,m3d_names,m3d_free_funcs,m3d_sum);
  ....
  M = m3d_get(3,4,5);
  ....
  /* making use of M->me[i][j][k], where i,j,k are non-negative and 
	i < 3, j < 4, k < 5 */
  ....
  mem_info_file(stdout,M3D_LIST);  /* info on the number of allocated 
				      bytes of memory for types 
	 			      on the list M3D_LIST */
  ....
  m3d_free(M);  /* if M is not necessary */
  ....
}


We can now use the function mem_info_file() for getting information about
the number of bytes of allocated memory and number of allocated variables
of type MAT3D; mem_stat_reg_list() for registering variables of this type
and mem_stat_mark() and mem_stat_free_list() for deallocating static
variables of this type.



In the similar way you can create you own list of errors and attach it to
the system. See the functions: 

  int err_list_attach(int list_num, int list_len, char **err_ptr,
		      int warn);  /* for attaching a list of errors */

  int err_is_list_attached(int list_num);  /* checking if a list 
                                                    is attached */

  extern  int err_list_free(int list_num);   /* freeing a list of errors */

where list_num is the number of the error list, list_len is the number of
errors on the list, err_ptr is the character string explaining the error
and warn can be TRUE if this is only a warning (the program continues to
run) or it can be FALSE if it is an error (the program stops).

The examples are the standard errors (error list 0) and warnings
(error list 1) which are in the file err.c


				David Stewart and Zbigniew Leyk, 1993