gmathlib.dox 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357
  1. /*! \page gmathlib GRASS Numerical math interface
  2. <!-- doxygenized from "GRASS 5 Programmer's Manual"
  3. by M. Neteler 8/2005
  4. -->
  5. <title>GRASS Numerical math interface</title>
  6. \section gmathintro Numerical math interface with optional support of LAPACK/BLAS
  7. <P>
  8. Author: David D. Gray<br>
  9. Additions: Brad Douglas
  10. <P>
  11. (under development)
  12. <BR>
  13. <P>
  14. This chapter provides an explanation of how numerical algebra routines from
  15. LAPACK/BLAS can be accessed through the GRASS GIS library "gmath". Most of
  16. the functions are wrapper modules for linear algebra problems, a few are locally
  17. implemented.
  18. <BR>
  19. <P>
  20. Getting BLAS/LAPACK (one package) if not already provided by the system:
  21. <BR><TT><A HREF="http://www.netlib.org/lapack/">http://www.netlib.org/lapack/</A></TT>
  22. <BR><TT><A HREF="http://netlib.bell-labs.com/netlib/master/readme.html">http://netlib.bell-labs.com/netlib/master/readme.html</A></TT>
  23. <BR>
  24. <P>
  25. Pre-compiled binaries of LAPACK/BLAS are provided on many Linux
  26. distributions.
  27. <P>
  28. \section Implementation Implementation
  29. <P>
  30. The function name convention is as follows:
  31. <P>
  32. <OL>
  33. <LI>G_matrix_*() : corresponding to Level3 BLAS (matrix -matrix) ,
  34. </LI>
  35. <LI>G_matvect_*() : corresponding to Level2 BLAS (vector-matrix) and
  36. </LI>
  37. <LI>G_vector_*() : corresponding to Level1 BLAS (vector-vector)
  38. </LI>
  39. </OL>
  40. <P>
  41. \section matrix_matrix_functions matrix -matrix functions
  42. <P>
  43. mat_struct *G_matrix_init (int rows, int cols, int ldim)<br>
  44. Initialise a matrix structure Initialise a matrix structure. Set
  45. the number of rows with
  46. the first parameter and columns with the second. The third parameter, lead
  47. dimension (>= no. of rows) needs attention by the programmer as it is
  48. related to the Fortran workspace:
  49. <P>
  50. A 3x3 matrix would be stored as
  51. <P>
  52. &nbsp;&nbsp;[ x x x _ ][ x x x _ ][ x x x _ ]
  53. <BR>
  54. <P>
  55. This work space corresponds to the sequence:
  56. <P>
  57. (1,1) (2,1) (3,1) and unused are (1,2) (2,2) ... and so on, ie. it is column major.
  58. So although the programmer uses the normal parameter sequence of (row, col)
  59. the entries traverse the data column by column instead of row by row. Each
  60. block in the workspace must be large enough (here 4) to hold a whole column
  61. (3) , and trailing spaces are unused. The number of rows (ie. size of a
  62. column) is therefore not necessarily the same size as the block size
  63. allocated to hold them (called the "lead dimension") . Some operations may
  64. produce matrices a different size from the inputs, but still write to the
  65. same workspace. This may seem awkward, but it does make for efficient code.
  66. Unfortunately, this creates a responsibility on the programmer to specify the
  67. lead dimension (>= no. of rows). In most cases the programmer can just use
  68. the rows. So for 3 rows/2 cols it would be called:
  69. <P>
  70. G_matrix_init (3, 2, 3);
  71. <P>
  72. mat_struct *G_matrix_set (mat_struct *A, int rows, int cols, int ldim)<br>
  73. Set parameters for a matrix structureSet parameters for a matrix
  74. structure that is allocated but not yet initialised fully.
  75. <P>
  76. <B>Note:</B>
  77. <BR>
  78. <P>
  79. G_matrix_set() is an alternative to G_matrix_init() . G_matrix_init()
  80. initialises and returns a pointer to a dynamically allocated matrix
  81. structure (ie. in the process heap) . G_matrix_set() sets the parameters for
  82. an already created matrix structure. In both cases the data workspace still
  83. has to be allocated dynamically.
  84. <P>
  85. Example 1:
  86. <P>
  87. \verbatim
  88. mat_struct *A;
  89. G_matrix_set (A, 4, 3, 4);
  90. \endverbatim
  91. <BR>
  92. <P>
  93. Example 2:
  94. <P>
  95. \verbatim
  96. mat_struct A; /* Allocated on the local stack */
  97. G_matrix_set (&A, 4, 3, 4) ;
  98. \endverbatim
  99. <BR>
  100. <P>
  101. mat_struct *G_matrix_add (mat_struct *mt1, mat_struct *mt2)<br>
  102. Add two matrices. Add two matrices and return the result. The receiving structure
  103. should not be initialised, as the matrix is created by the routine.
  104. <P>
  105. mat_struct *G_matrix_product (mat_struct *mt1, mat_struct *mt2)<br>
  106. Multiply two matricesMultiply two matrices and return the result.
  107. The receiving structure should not be initialised, as the matrix is created
  108. by the routine.
  109. <P>
  110. mat_struct *G_matrix_scale (mat_struct *mt1, const double c)<br>
  111. Scale matrix Scale the matrix by a given scalar value and return the
  112. result. The receiving structure should not be initialised, as the matrix is
  113. created by the routine.
  114. <P>
  115. mat_struct *G_matrix_subtract (mat_struct *mt1, mat_struct *mt2)<br>
  116. Subtract two matrices. Subtract two matrices and return the result.
  117. The receiving structure should not be initialised, as the matrix is created
  118. by the routine.
  119. <P>
  120. mat_struct *G_matrix_copy (const mat_struct *A)<br>
  121. Copy a matrix. Copy a matrix by exactly duplicating its structure.
  122. <P>
  123. mat_struct *G_matrix_transpose (mat_struct *mt)<br>
  124. Transpose a matrix. Transpose a matrix by creating a new one and
  125. populating with transposed element s.
  126. <P>
  127. void G_matrix_print (mat_struct *mt)<br>
  128. Print out a matrix. Print out a representation of the matrix to
  129. standard output.
  130. <P>
  131. int G_matrix_LU_solve (const mat_struct *mt1, mat_struct **xmat0, const
  132. mat_struct *bmat, mat_type mtype)<br>
  133. Solve a general system A.X=B. Solve a general
  134. system A.X=B, where A is a NxN matrix, X and B are NxC matrices, and we are to
  135. solve for C arrays in X given B. Uses LU decomposition.
  136. <BR>
  137. <P>
  138. Links to LAPACK function dgesv_() and similar to perform the core routine.
  139. (By default solves for a general non-symmetric matrix .)
  140. <P>
  141. mtype is a flag to indicate what kind of matrix (real/complex, Hermitian,
  142. symmetric, general etc.) is used (NONSYM, SYM, HERMITIAN) .
  143. <P>
  144. <B>Warning:</B> NOT YET COMPLETE: only some solutions' options
  145. available. Now, only general real matrix is supported.
  146. <P>
  147. mat_struct *G_matrix_inverse (mat_struct *mt)<br>
  148. Matrix inversion using
  149. LU decomposition Calls G_matrix_LU_solve() to obtain matrix inverse using
  150. LU decomposition. Returns NULL on failure.
  151. <P>
  152. void G_matrix_free (mat_struct *mt)<br> Free up allocated matrix Free up
  153. allocated matrix.
  154. <P>
  155. int G_matrix_set_element (mat_struct *mt, int rowval, int colval, double val)<br>
  156. Set the value of the (i,j) th element Set the value of the
  157. (i,j) th element to a double value. Index values are C-like ie. zero-based.
  158. The row number is given first as is conventional. Returns -1 if the
  159. accessed cell is outside the bounds.
  160. <P>
  161. double G_matrix_get_element (mat_struct *mt, int rowval, int colval)<br>
  162. Retrieve value of the (i,j) th element Retrieve the value of the
  163. (i,j) th element to a double value. Index values are C-like ie. zero-based.
  164. <P>
  165. <B>Note:</B> Does currently not set an error flag for bounds checking.
  166. <P>
  167. \section matrix_Vector_functions Matrix-Vector functions
  168. <P>
  169. vec_struct *G_matvect_get_column (mat_struct *mt, int
  170. col) Retrieve a column of matrix Retrieve a column of the matrix to a vector
  171. structure. The receiving structure should not be initialised, as the
  172. vector is created by the routine. Col 0 is the first column.
  173. <P>
  174. vec_struct *G_matvect_get_row (mat_struct *mt, int
  175. row) Retrieve a row of matrix Retrieve a row of the matrix to a vector
  176. structure. The receiving structure should not be initialised, as the
  177. vector is created by the routine. Row 0 is the first number.
  178. <P>
  179. int G_matvect_extract_vector (mat_struct *mt, vtype vt, int
  180. indx) Convert matrix to vectorConvert the current matrix structure to
  181. a vector structure. The vtype is RVEC or CVEC which specifies a row vector or
  182. column vector. The indx indicates the row/column number (zero based) .
  183. <P>
  184. int G_matvect_retrieve_matrix (vec_struct *vc) Revert a
  185. vector to matrix Revert a vector structure to a matrix .
  186. <P>
  187. \section Vector_Vector_functions Vector-Vector functions
  188. vec_struct *G_vector_init (int cells, int ldim, vtype vt) Initialise
  189. a vector structure Initialise a vector structure. The vtype is RVEC or
  190. CVEC which specifies a row vector or column vector.
  191. <P>
  192. int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int vindx) Set
  193. parameters for vector structureSet parameters for a vector structure that is
  194. allocated but not yet initialised fully. The vtype is RVEC or
  195. CVEC which specifies a row vector or column vector.
  196. <P>
  197. vec_struct *G_vector_copy (const vec_struct *vc1, int
  198. comp_flag) Copy a vectorCopy a vector to a new vector structure. This preserves
  199. the underlying structure unless you specify DO_COMPACT comp_flag:
  200. <P>
  201. 0 Eliminate unnecessary rows (cols) in matrix
  202. <BR>
  203. 1 ... or not
  204. <P>
  205. double G_vector_norm_euclid (vec_struct *vc) Calculates euclidean
  206. norm Calculates the euclidean norm of a row or column vector, using BLAS
  207. routine dnrm2_()
  208. <P>
  209. double G_vector_norm_maxval (vec_struct *vc, int vflag) Calculates
  210. maximum value Calculates the maximum value of a row or column vector.
  211. The vflag setting defines which value to be calculated:
  212. <P>
  213. vflag:
  214. <P>
  215. 1 Indicates maximum value
  216. <BR> -1 Indicates minimum value
  217. <BR>
  218. 0 Indicates absolute value [???]
  219. <P>
  220. <B>Note:</B> More functions and wrappers will be implemented soon (11/2000) .
  221. <P>
  222. \section Notes Notes
  223. The numbers of rows and columns is stored in the matrix structure:
  224. <P>
  225. \verbatim
  226. G_message (" M1 rows: %d, M1 cols: %d", m1->rows, m1->cols);
  227. \endverbatim
  228. <P>
  229. Draft Notes:
  230. <P>
  231. * G_vector_free() can be done by G_matrix_free() .
  232. <P>
  233. \verbatim
  234. #define G_vector_free(x) G_matrix_free( (x) )
  235. \endverbatim
  236. <P>
  237. * Ax calculations can be done with G_matrix_multiply()
  238. <P>
  239. * Vector print can be done by G_matrix_print() .
  240. <P>
  241. \verbatim
  242. #define G_vector_print(x) G_matrix_print( (x) )
  243. \endverbatim
  244. <P>
  245. \section Example Example
  246. The Makefile needs a reference to $(GMATHLIB) in LIBES line.
  247. <P>
  248. Example module:
  249. <P>
  250. \verbatim
  251. #include <grass/config.h>
  252. #include <grass/gis.h>
  253. #include <grass/gmath.h>
  254. int
  255. main (int argc, char *argv[])
  256. {
  257. mat_struct *matrix;
  258. double val;
  259. /* init a 3x2 matrix */
  260. matrix = G_matrix_init (3, 2, 3);
  261. /* set cell element 0,0 in matrix to 2.2: */
  262. G_matrix_set_element (matrix, 0, 0, 2.2);
  263. /* retrieve this value */
  264. val = G_matrix_get_element (matrix, 0, 0);
  265. /* print it */
  266. G_message ( "Value: %g", val);
  267. /* Free the memory */
  268. G_matrix_free (matrix);
  269. return 0;
  270. }
  271. \endverbatim
  272. <P>
  273. See \ref Compiling_and_Installing_GRASS_Modules for a complete
  274. discussion of Makefiles.
  275. */