|
GRASS Programmer's Manual
6.5.svn(2014)-r66266
|
|
GRASS Programmer's Manual
6.5.svn(2014)-r66266
|
Author: David D. Gray
Additions: Brad Douglas
(under development)
This chapter provides an explanation of how numerical algebra routines from LAPACK/BLAS can be accessed through the GRASS GIS library "gmath". Most of the functions are wrapper modules for linear algebra problems, a few are locally implemented.
Getting BLAS/LAPACK (one package) if not already provided by the system:
http://www.netlib.org/lapack/
http://netlib.bell-labs.com/netlib/master/readme.html
Pre-compiled binaries of LAPACK/BLAS are provided on many Linux distributions.
The function name convention is as follows:
mat_struct *G_matrix_init (int rows, int cols, int ldim)
Initialise a matrix structure Initialise a matrix structure. Set the number of rows with the first parameter and columns with the second. The third parameter, lead dimension (>= no. of rows) needs attention by the programmer as it is related to the Fortran workspace:
A 3x3 matrix would be stored as
[ x x x _ ][ x x x _ ][ x x x _ ]
This work space corresponds to the sequence:
(1,1) (2,1) (3,1) and unused are (1,2) (2,2) ... and so on, ie. it is column major. So although the programmer uses the normal parameter sequence of (row, col) the entries traverse the data column by column instead of row by row. Each block in the workspace must be large enough (here 4) to hold a whole column (3) , and trailing spaces are unused. The number of rows (ie. size of a column) is therefore not necessarily the same size as the block size allocated to hold them (called the "lead dimension") . Some operations may produce matrices a different size from the inputs, but still write to the same workspace. This may seem awkward, but it does make for efficient code. Unfortunately, this creates a responsibility on the programmer to specify the lead dimension (>= no. of rows). In most cases the programmer can just use the rows. So for 3 rows/2 cols it would be called:
G_matrix_init (3, 2, 3);
mat_struct *G_matrix_set (mat_struct *A, int rows, int cols, int ldim)
Set parameters for a matrix structureSet parameters for a matrix structure that is allocated but not yet initialised fully.
Note:
G_matrix_set() is an alternative to G_matrix_init() . G_matrix_init() initialises and returns a pointer to a dynamically allocated matrix structure (ie. in the process heap) . G_matrix_set() sets the parameters for an already created matrix structure. In both cases the data workspace still has to be allocated dynamically.
Example 1:
mat_struct *A; G_matrix_set (A, 4, 3, 4);
Example 2:
mat_struct A; /* Allocated on the local stack */ G_matrix_set (&A, 4, 3, 4) ;
mat_struct *G_matrix_add (mat_struct *mt1, mat_struct *mt2)
Add two matrices. Add two matrices and return the result. The receiving structure should not be initialised, as the matrix is created by the routine.
mat_struct *G_matrix_product (mat_struct *mt1, mat_struct *mt2)
Multiply two matricesMultiply two matrices and return the result. The receiving structure should not be initialised, as the matrix is created by the routine.
mat_struct *G_matrix_scale (mat_struct *mt1, const double c)
Scale matrix Scale the matrix by a given scalar value and return the result. The receiving structure should not be initialised, as the matrix is created by the routine.
mat_struct *G_matrix_subtract (mat_struct *mt1, mat_struct *mt2)
Subtract two matrices. Subtract two matrices and return the result. The receiving structure should not be initialised, as the matrix is created by the routine.
mat_struct *G_matrix_copy (const mat_struct *A)
Copy a matrix. Copy a matrix by exactly duplicating its structure.
mat_struct *G_matrix_transpose (mat_struct *mt)
Transpose a matrix. Transpose a matrix by creating a new one and populating with transposed element s.
void G_matrix_print (mat_struct *mt)
Print out a matrix. Print out a representation of the matrix to standard output.
int G_matrix_LU_solve (const mat_struct *mt1, mat_struct **xmat0, const mat_struct *bmat, mat_type mtype)
Solve a general system A.X=B. Solve a general system A.X=B, where A is a NxN matrix, X and B are NxC matrices, and we are to solve for C arrays in X given B. Uses LU decomposition.
Links to LAPACK function dgesv_() and similar to perform the core routine. (By default solves for a general non-symmetric matrix .)
mtype is a flag to indicate what kind of matrix (real/complex, Hermitian, symmetric, general etc.) is used (NONSYM, SYM, HERMITIAN) .
Warning: NOT YET COMPLETE: only some solutions' options available. Now, only general real matrix is supported.
mat_struct *G_matrix_inverse (mat_struct *mt)
Matrix inversion using LU decomposition Calls G_matrix_LU_solve() to obtain matrix inverse using LU decomposition. Returns NULL on failure.
void G_matrix_free (mat_struct *mt)
Free up allocated matrix Free up allocated matrix.
int G_matrix_set_element (mat_struct *mt, int rowval, int colval, double val)
Set the value of the (i,j) th element Set the value of the (i,j) th element to a double value. Index values are C-like ie. zero-based. The row number is given first as is conventional. Returns -1 if the accessed cell is outside the bounds.
double G_matrix_get_element (mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j) th element Retrieve the value of the (i,j) th element to a double value. Index values are C-like ie. zero-based.
Note: Does currently not set an error flag for bounds checking.
vec_struct *G_matvect_get_column (mat_struct *mt, int col) Retrieve a column of matrix Retrieve a column of the matrix to a vector structure. The receiving structure should not be initialised, as the vector is created by the routine. Col 0 is the first column.
vec_struct *G_matvect_get_row (mat_struct *mt, int row) Retrieve a row of matrix Retrieve a row of the matrix to a vector structure. The receiving structure should not be initialised, as the vector is created by the routine. Row 0 is the first number.
int G_matvect_extract_vector (mat_struct *mt, vtype vt, int indx) Convert matrix to vectorConvert the current matrix structure to a vector structure. The vtype is RVEC or CVEC which specifies a row vector or column vector. The indx indicates the row/column number (zero based) .
int G_matvect_retrieve_matrix (vec_struct *vc) Revert a vector to matrix Revert a vector structure to a matrix .
vec_struct *G_vector_init (int cells, int ldim, vtype vt) Initialise a vector structure Initialise a vector structure. The vtype is RVEC or CVEC which specifies a row vector or column vector.
int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int vindx) Set parameters for vector structureSet parameters for a vector structure that is allocated but not yet initialised fully. The vtype is RVEC or CVEC which specifies a row vector or column vector.
vec_struct *G_vector_copy (const vec_struct *vc1, int comp_flag) Copy a vectorCopy a vector to a new vector structure. This preserves the underlying structure unless you specify DO_COMPACT comp_flag:
0 Eliminate unnecessary rows (cols) in matrix
1 ... or not
double G_vector_norm_euclid (vec_struct *vc) Calculates euclidean norm Calculates the euclidean norm of a row or column vector, using BLAS routine dnrm2_()
double G_vector_norm_maxval (vec_struct *vc, int vflag) Calculates maximum value Calculates the maximum value of a row or column vector. The vflag setting defines which value to be calculated:
vflag:
1 Indicates maximum value
-1 Indicates minimum value
0 Indicates absolute value [???]
Note: More functions and wrappers will be implemented soon (11/2000) .
The numbers of rows and columns is stored in the matrix structure:
G_message (" M1 rows: %d, M1 cols: %d", m1->rows, m1->cols);
Draft Notes:
G_vector_free() can be done by G_matrix_free() .
#define G_vector_free(x) G_matrix_free( (x) )
Ax calculations can be done with G_matrix_multiply()
Vector print can be done by G_matrix_print() .
#define G_vector_print(x) G_matrix_print( (x) )
The Makefile needs a reference to in LIBES line.
Example module:
#include <grass/config.h>
#include <grass/gis.h>
#include <grass/gmath.h>
int
main (int argc, char *argv[])
{
mat_struct *matrix;
double val;
/* init a 3x2 matrix */
matrix = G_matrix_init (3, 2, 3);
/* set cell element 0,0 in matrix to 2.2: */
G_matrix_set_element (matrix, 0, 0, 2.2);
/* retrieve this value */
val = G_matrix_get_element (matrix, 0, 0);
/* print it */
G_message ( "Value: %g", val);
/* Free the memory */
G_matrix_free (matrix);
return 0;
}
See Compiling and Installing GRASS Modules for a complete discussion of Makefiles.
1.8.5