GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
ccmath_grass_wrapper.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass numerical math interface
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> googlemail <dot> com
7 *
8 * PURPOSE: ccmath library function wrapper
9 * part of the gmath library
10 *
11 * COPYRIGHT: (C) 2010 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 
19 
20 #if defined(HAVE_CCMATH)
21 #include <ccmath.h>
22 #else
23 #include <grass/ccmath_grass.h>
24 #endif
25 /**
26  * Documentation and ccmath library version 2.2.1 by Daniel A. Atkinson
27  *
28  Chapter 1
29 
30  LINEAR ALGEBRA
31 
32  Summary
33 
34  The matrix algebra library contains functions that
35  perform the standard computations of linear algebra.
36  General areas covered are:
37 
38  o Solution of Linear Systems
39  o Matrix Inversion
40  o Eigensystem Analysis
41  o Matrix Utility Operations
42  o Singular Value Decomposition
43 
44  The operations covered here are fundamental to many
45  areas of mathematics and statistics. Thus, functions
46  in this library segment are called by other library
47  functions. Both real and complex valued matrices
48  are covered by functions in the first four of these
49  categories.
50 
51 
52  Notes on Contents
53 
54  Functions in this library segment provide the basic operations of
55  numerical linear algebra and some useful utility functions for operations on
56  vectors and matrices. The following list describes the functions available for
57  operations with real-valued matrices.
58 
59 
60  o Solving and Inverting Linear Systems:
61 
62  solv --------- solve a general system of real linear equations.
63  solvps ------- solve a real symmetric linear system.
64  solvru ------- solve a real right upper triangular linear system.
65  solvtd ------- solve a tridiagonal real linear system.
66 
67  minv --------- invert a general real square matrix.
68  psinv -------- invert a real symmetric matrix.
69  ruinv -------- invert a right upper triangular matrix.
70 
71 
72  The solution of a general linear system and efficient algorithms for
73  solving special systems with symmetric and tridiagonal matrices are provided
74  by these functions. The general solution function employs a LU factorization
75  with partial pivoting and it is very robust. It will work efficiently on any
76  problem that is not ill-conditioned. The symmetric matrix solution is based
77  on a modified Cholesky factorization. It is best used on positive definite
78  matrices that do not require pivoting for numeric stability. Tridiagonal
79  solvers require order-N operations (N = dimension). Thus, they are highly
80  recommended for this important class of sparse systems. Two matrix inversion
81  routines are provided. The general inversion function is again LU based. It
82  is suitable for use on any stable (ie. well-conditioned) problem. The
83  Cholesky based symmetric matrix inversion is efficient and safe for use on
84  matrices known to be positive definite, such as the variance matrices
85  encountered in statistical computations. Both the solver and the inverse
86  functions are designed to enhance data locality. They are very effective
87  on modern microprocessors.
88 
89 
90  o Eigensystem Analysis:
91 
92  eigen ------ extract all eigen values and vectors of a real
93  symmetric matrix.
94  eigval ----- extract the eigen values of a real symmetric matrix.
95  evmax ------ compute the eigen value of maximum absolute magnitude
96  and its corresponding vector for a symmetric matrix.
97 
98 
99  Eigensystem functions operate on real symmetric matrices. Two forms of
100  the general eigen routine are provided because the computation of eigen values
101  only is much faster when vectors are not required. The basic algorithms use
102  a Householder reduction to tridiagonal form followed by QR iterations with
103  shifts to enhance convergence. This has become the accepted standard for
104  symmetric eigensystem computation. The evmax function uses an efficient
105  iterative power method algorithm to extract the eigen value of maximum
106  absolute size and the corresponding eigenvector.
107 
108 
109  o Singular Value Decomposition:
110 
111  svdval ----- compute the singular values of a m by n real matrix.
112  sv2val ----- compute the singular values of a real matrix
113  efficiently for m >> n.
114  svduv ------ compute the singular values and the transformation
115  matrices u and v for a real m by n matrix.
116  sv2uv ------ compute the singular values and transformation
117  matrices efficiently for m >> n.
118  svdu1v ----- compute the singular values and transformation
119  matrices u1 and v, where u1 overloads the input
120  with the first n column vectors of u.
121  sv2u1v ----- compute the singular values and the transformation
122  matrices u1 and v efficiently for m >> n.
123 
124 
125  Singular value decomposition is extremely useful when dealing with linear
126  systems that may be singular. Singular values with values near zero are flags
127  of a potential rank deficiency in the system matrix. They can be used to
128  identify the presence of an ill-conditioned problem and, in some cases, to
129  deal with the potential instability. They are applied to the linear least
130  squares problem in this library. Singular values also define some important
131  matrix norm parameters such as the 2-norm and the condition value. A complete
132  decomposition provides both singular values and an orthogonal decomposition of
133  vector spaces related to the matrix identifying the range and null-space.
134  Fortunately, a highly stable algorithm based on Householder reduction to
135  bidiagonal form and QR rotations can be used to implement the decomposition.
136  The library provides two forms with one more efficient when the dimensions
137  satisfy m > (3/2)n.
138 
139  General Technical Comments
140 
141  Efficient computation with matrices on modern processors must be
142  adapted to the storage scheme employed for matrix elements. The functions
143  of this library segment do not employ the multidimensional array intrinsic
144  of the C language. Access to elements employs the simple row-major scheme
145  described here.
146 
147  Matrices are modeled by the library functions as arrays with elements
148  stored in row order. Thus, the element in the jth row and kth column of
149  the n by n matrix M, stored in the array mat[], is addressed by
150 
151  M[j,k] = mat[n*j+k] , with 0 =< j,k <= n-1 .
152 
153  (Remember that C employs zero as the starting index.) The storage order has
154  important implications for data locality.
155 
156  The algorithms employed here all have excellent numerical stability, and
157  the default double precision arithmetic of C enhances this. Thus, any
158  problems encountered in using the matrix algebra functions will almost
159  certainly be due to an ill-conditioned matrix. (The Hilbert matrices,
160 
161  H[i,j] = 1/(1+i+j) for i,j < n
162 
163  form a good example of such ill-conditioned systems.) We remind the reader
164  that the appropriate response to such ill-conditioning is to seek an
165  alternative approach to the problem. The option of increasing precision has
166  already been exploited. Modification of the linear algebra algorithm code is
167  not normally effective in an ill-conditioned problem.
168 
169 ------------------------------------------------------------------------------
170  FUNCTION SYNOPSES
171 ------------------------------------------------------------------------------
172 
173  Linear System Solutions:
174 -----------------------------------------------------------------------------
175 */
176 /**
177  \brief Solve a general linear system A*x = b.
178 
179  \param a = array containing system matrix A in row order (altered to L-U factored form by computation)
180  \param b = array containing system vector b at entry and solution vector x at exit
181  \param n = dimension of system
182  \return 0 -> normal exit; -1 -> singular input
183  */
184 int G_math_solv(double **a,double *b,int n)
185 {
186  return solv(a[0],b, n);
187 }
188 
189 
190 /**
191  \brief Solve a symmetric positive definite linear system S*x = b.
192 
193  \param a = array containing system matrix S (altered to Cholesky upper right factor by computation)
194  \param b = array containing system vector b as input and solution vector x as output
195  \param n = dimension of system
196  \return: 0 -> normal exit; -1 -> input matrix not positive definite
197  */
198  int G_math_solvps(double **a,double *b,int n)
199 {
200  return solvps(a[0], b,n);
201 }
202 
203 
204 /**
205  \brief Solve a tridiagonal linear system M*x = y.
206 
207  \param a = array containing m+1 diagonal elements of M
208  \param b = array of m elements below the main diagonal of M
209  \param c = array of m elements above the main diagonal
210  \param x = array containing the system vector y initially, and the solution vector at exit (m+1 elements)
211  \param m = dimension parameter ( M is (m+1)x(m+1) )
212 
213 */
214 void G_math_solvtd(double *a,double *b,double *c,double *x,int m)
215 {
216  solvtd(a, b, c, x, m);
217  return;
218 }
219 
220 
221 /*
222  \brief Solve an upper right triangular linear system T*x = b.
223 
224  \param a = pointer to array of upper right triangular matrix T
225  \param b = pointer to array of system vector The computation overloads this with the solution vector x.
226  \param n = dimension (dim(a)=n*n,dim(b)=n)
227  \return value: f = status flag, with 0 -> normal exit, -1 -> system singular
228 */
229 int G_math_solvru(double **a,double *b,int n)
230 {
231  return solvru(a[0], b, n);
232 }
233 
234 
235 /**
236  \brief Invert (in place) a general real matrix A -> Inv(A).
237 
238  \param a = array containing the input matrix A. This is converted to the inverse matrix.
239  \param n = dimension of the system (i.e. A is n x n )
240  \return: 0 -> normal exit, 1 -> singular input matrix
241 */
242 int G_math_minv(double **a,int n)
243 {
244  return minv(a[0], n);
245 }
246 
247 
248 /**
249  \brief Invert (in place) a symmetric real matrix, V -> Inv(V).
250 
251  The input matrix V is symmetric (V[i,j] = V[j,i]).
252  \param a = array containing a symmetric input matrix. This is converted to the inverse matrix.
253  \param n = dimension of the system (dim(v)=n*n)
254  \return: 0 -> normal exit 1 -> input matrix not positive definite
255 */
256 int G_math_psinv(double **a,int n)
257 {
258  return psinv( a[0], n);
259 }
260 
261 
262 /**
263  \brief Invert an upper right triangular matrix T -> Inv(T).
264 
265  \param a = pointer to array of upper right triangular matrix, This is replaced by the inverse matrix.
266  \param n = dimension (dim(a)=n*n)
267  \return value: status flag, with 0 -> matrix inverted -1 -> matrix singular
268 */
269 int G_math_ruinv(double **a,int n)
270 {
271  return ruinv(a[0], n);
272 }
273 
274 
275 /*
276 -----------------------------------------------------------------------------
277 
278  Symmetric Eigensystem Analysis:
279 -----------------------------------------------------------------------------
280 */
281 /**
282 
283  \brief Compute the eigenvalues of a real symmetric matrix A.
284 
285  \param a = pointer to array of symmetric n by n input matrix A. The computation alters these values.
286  \param ev = pointer to array of the output eigenvalues
287  \param n = dimension parameter (dim(a)= n*n, dim(ev)= n)
288 */
289 void G_math_eigval(double **a,double *ev,int n)
290 {
291  eigval(a[0], ev, n);
292  return;
293 }
294 
295 
296 /**
297  \brief Compute the eigenvalues and eigenvectors of a real symmetric matrix A.
298 
299  The input and output matrices are related by
300 
301  A = E*D*E~ where D is the diagonal matrix of eigenvalues
302  D[i,j] = ev[i] if i=j and 0 otherwise.
303 
304  The columns of E are the eigenvectors.
305 
306  \param a = pointer to store for symmetric n by n input matrix A. The computation overloads this with an orthogonal matrix of eigenvectors E.
307  \param ev = pointer to the array of the output eigenvalues
308  \param n = dimension parameter (dim(a)= n*n, dim(ev)= n)
309 */
310 void G_math_eigen(double **a,double *ev,int n)
311 {
312  eigen(a[0], ev, n);
313  return;
314 }
315 
316 
317 /*
318  \brief Compute the maximum (absolute) eigenvalue and corresponding eigenvector of a real symmetric matrix A.
319 
320 
321  \param a = array containing symmetric input matrix A
322  \param u = array containing the n components of the eigenvector at exit (vector normalized to 1)
323  \param n = dimension of system
324  \return: ev = eigenvalue of A with maximum absolute value HUGE -> convergence failure
325 */
326 double G_math_evmax(double **a,double *u,int n)
327 {
328  return evmax(a[0], u, n);
329 }
330 
331 
332 /*
333 ------------------------------------------------------------------------------
334 
335  Singular Value Decomposition:
336 ------------------------------------------------------------------------------
337 
338  A number of versions of the Singular Value Decomposition (SVD)
339  are implemented in the library. They support the efficient
340  computation of this important factorization for a real m by n
341  matrix A. The general form of the SVD is
342 
343  A = U*S*V~ with S = | D |
344  | 0 |
345 
346  where U is an m by m orthogonal matrix, V is an n by n orthogonal matrix,
347  D is the n by n diagonal matrix of singular value, and S is the singular
348  m by n matrix produced by the transformation.
349 
350  The singular values computed by these functions provide important
351  information on the rank of the matrix A, and on several matrix
352  norms of A. The number of non-zero singular values d[i] in D
353  equal to the rank of A. The two norm of A is
354 
355  ||A|| = max(d[i]) , and the condition number is
356 
357  k(A) = max(d[i])/min(d[i]) .
358 
359  The Frobenius norm of the matrix A is
360 
361  Fn(A) = Sum(i=0 to n-1) d[i]^2 .
362 
363  Singular values consistent with zero are easily recognized, since
364  the decomposition algorithms have excellent numerical stability.
365  The value of a 'zero' d[i] is no larger than a few times the
366  computational rounding error e.
367 
368  The matrix U1 is formed from the first n orthonormal column vectors
369  of U. U1[i,j] = U[i,j] for i = 1 to m and j = 1 to n. A singular
370  value decomposition of A can also be expressed in terms of the m by\
371  n matrix U1, with
372 
373  A = U1*D*V~ .
374 
375  SVD functions with three forms of output are provided. The first
376  form computes only the singular values, while the second computes
377  the singular values and the U and V orthogonal transformation
378  matrices. The third form of output computes singular values, the
379  V matrix, and saves space by overloading the input array with
380  the U1 matrix.
381 
382  Two forms of decomposition algorithm are available for each of the
383  three output types. One is computationally efficient when m ~ n.
384  The second, distinguished by the prefix 'sv2' in the function name,
385  employs a two stage Householder reduction to accelerate computation
386  when m substantially exceeds n. Use of functions of the second form
387  is recommended for m > 2n.
388 
389  Singular value output from each of the six SVD functions satisfies
390 
391  d[i] >= 0 for i = 0 to n-1.
392 -------------------------------------------------------------------------------
393 */
394 
395 
396 /**
397  \brief Compute the singular values of a real m by n matrix A.
398 
399 
400  \param d = pointer to double array of dimension n (output = singular values of A)
401  \param a = pointer to store of the m by n input matrix A (A is altered by the computation)
402  \param m = number of rows in A
403  \param n = number of columns in A (m>=n required)
404  \return value: status flag with: 0 -> success -1 -> input error m < n
405 
406 */
407 int G_math_svdval(double *d,double **a,int m,int n)
408 {
409  return svdval(d, a[0], m, n);
410 }
411 
412 
413 /**
414 
415  \brief Compute singular values when m >> n.
416 
417  \param d = pointer to double array of dimension n (output = singular values of A)
418  \param a = pointer to store of the m by n input matrix A (A is altered by the computation)
419  \param m = number of rows in A
420  \param n = number of columns in A (m>=n required)
421  \return value: status flag with: 0 -> success -1 -> input error m < n
422 */
423 int G_math_sv2val(double *d,double **a,int m,int n)
424 {
425  return sv2val(d, a[0], m, n);
426 }
427 
428 
429 /*
430  \brief Compute the singular value transformation S = U~*A*V.
431 
432  \param d = pointer to double array of dimension n (output = singular values of A)
433  \param a = pointer to store of the m by n input matrix A (A is altered by the computation)
434  \param u = pointer to store for m by m orthogonal matrix U
435  \param v = pointer to store for n by n orthogonal matrix V
436  \param m = number of rows in A
437  \param n = number of columns in A (m>=n required)
438  \return value: status flag with: 0 -> success -1 -> input error m < n
439 */
440 int G_math_svduv(double *d,double **a,double **u,int m,double **v,int n)
441 {
442  return svduv(d, a[0], u[0], m, v[0], n);
443 }
444 
445 
446 /**
447  \brief Compute the singular value transformation when m >> n.
448 
449  \param d = pointer to double array of dimension n (output = singular values of A)
450  \param a = pointer to store of the m by n input matrix A (A is altered by the computation)
451  \param u = pointer to store for m by m orthogonal matrix U
452  \param v = pointer to store for n by n orthogonal matrix V
453  \param m = number of rows in A
454  \param n = number of columns in A (m>=n required)
455  \return value: status flag with: 0 -> success -1 -> input error m < n
456 */
457 int G_math_sv2uv(double *d,double **a,double **u,int m,double **v,int n)
458 {
459  return sv2uv(d, a[0], u[0], m, v[0], n);
460 }
461 
462 
463 /**
464 
465  \brief Compute the singular value transformation with A overloaded by the partial U-matrix.
466 
467  \param d = pointer to double array of dimension n
468  (output = singular values of A)
469  \param a = pointer to store of the m by n input matrix A (At output a is overloaded by the matrix U1 whose n columns are orthogonal vectors equal to the first n columns of U.)
470  \param v = pointer to store for n by n orthogonal matrix V
471  \param m = number of rows in A
472  \param n = number of columns in A (m>=n required)
473  \return value: status flag with: 0 -> success -1 -> input error m < n
474 
475 */
476 int G_math_svdu1v(double *d,double **a,int m,double **v,int n)
477 {
478  return svdu1v(d, a[0], m, v[0], n);
479 }
int svduv(double *d, double *a, double *u, int m, double *v, int n)
Definition: svduv.c:10
int G_math_solvru(double **a, double *b, int n)
int solvps(double *s, double *x, int n)
Definition: solvps.c:9
int G_math_solvps(double **a, double *b, int n)
Solve a symmetric positive definite linear system S*x = b.
void G_math_eigen(double **a, double *ev, int n)
Compute the eigenvalues and eigenvectors of a real symmetric matrix A.
int svdval(double *d, double *a, int m, int n)
Definition: svdval.c:10
double G_math_evmax(double **a, double *u, int n)
void eigen(double *a, double *eval, int n)
Definition: eigen.c:10
void G_math_eigval(double **a, double *ev, int n)
Compute the eigenvalues of a real symmetric matrix A.
int G_math_sv2uv(double *d, double **a, double **u, int m, double **v, int n)
Compute the singular value transformation when m >> n.
void solvtd(double *a, double *b, double *c, double *x, int m)
Definition: solvtd.c:8
int G_math_sv2val(double *d, double **a, int m, int n)
Compute singular values when m >> n.
#define x
void G_math_solvtd(double *a, double *b, double *c, double *x, int m)
Solve a tridiagonal linear system M*x = y.
int sv2val(double *d, double *a, int m, int n)
Definition: sv2val.c:10
int minv(double *a, int n)
Definition: minv.c:10
double b
Definition: r_raster.c:39
int psinv(double *v, int n)
Definition: psinv.c:9
int G_math_svduv(double *d, double **a, double **u, int m, double **v, int n)
int G_math_psinv(double **a, int n)
Invert (in place) a symmetric real matrix, V -> Inv(V).
int sv2uv(double *d, double *a, double *u, int m, double *v, int n)
Definition: sv2uv.c:10
int solvru(double *a, double *b, int n)
Definition: solvru.c:8
int svdu1v(double *d, double *a, int m, double *v, int n)
Definition: svdu1v.c:10
double evmax(double *a, double *u, int n)
Definition: evmax.c:10
int G_math_svdval(double *d, double **a, int m, int n)
Compute the singular values of a real m by n matrix A.
int G_math_svdu1v(double *d, double **a, int m, double **v, int n)
Compute the singular value transformation with A overloaded by the partial U-matrix.
int G_math_ruinv(double **a, int n)
Invert an upper right triangular matrix T -> Inv(T).
int G_math_minv(double **a, int n)
Invert (in place) a general real matrix A -> Inv(A).
int ruinv(double *a, int n)
Definition: ruinv.c:8
int solv(double *a, double *b, int n)
Definition: solv.c:10
int G_math_solv(double **a, double *b, int n)
Solve a general linear system A*x = b.
void eigval(double *a, double *eval, int n)
Definition: eigval.c:10