GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 */
184int 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*/
214void 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*/
229int 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*/
242int 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*/
256int 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*/
269int 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*/
289void 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*/
310void 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*/
326double 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*/
407int 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*/
423int 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*/
440int 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*/
457int 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*/
476int 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}
void solvtd(double *a, double *b, double *c, double *x, int m)
Definition: solvtd.c:8
int svdu1v(double *d, double *a, int m, double *v, int n)
Definition: svdu1v.c:10
void eigval(double *a, double *eval, int n)
Definition: eigval.c:10
int solv(double *a, double *b, int n)
Definition: solv.c:10
int sv2uv(double *d, double *a, double *u, int m, double *v, int n)
Definition: sv2uv.c:10
int svduv(double *d, double *a, double *u, int m, double *v, int n)
Definition: svduv.c:10
int psinv(double *v, int n)
Definition: psinv.c:9
int solvps(double *s, double *x, int n)
Definition: solvps.c:9
int ruinv(double *a, int n)
Definition: ruinv.c:8
double evmax(double *a, double *u, int n)
Definition: evmax.c:10
void eigen(double *a, double *eval, int n)
Definition: eigen.c:10
int svdval(double *d, double *a, int m, int n)
Definition: svdval.c:10
int solvru(double *a, double *b, int n)
Definition: solvru.c:8
int minv(double *a, int n)
Definition: minv.c:10
int sv2val(double *d, double *a, int m, int n)
Definition: sv2val.c:10
int G_math_solvru(double **a, double *b, int n)
int G_math_psinv(double **a, int n)
Invert (in place) a symmetric real matrix, V -> Inv(V).
int G_math_solv(double **a, double *b, int n)
Solve a general linear system A*x = b.
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_minv(double **a, int n)
Invert (in place) a general real matrix A -> Inv(A).
void G_math_eigen(double **a, double *ev, int n)
Compute the eigenvalues and eigenvectors of a real symmetric matrix A.
int G_math_svduv(double *d, double **a, double **u, int m, double **v, int n)
void G_math_eigval(double **a, double *ev, int n)
Compute the eigenvalues of a real symmetric matrix A.
int G_math_solvps(double **a, double *b, int n)
Solve a symmetric positive definite linear system S*x = b.
double G_math_evmax(double **a, double *u, int n)
void G_math_solvtd(double *a, double *b, double *c, double *x, int m)
Solve a tridiagonal linear system M*x = y.
int G_math_ruinv(double **a, int n)
Invert an upper right triangular matrix T -> Inv(T).
int G_math_sv2val(double *d, double **a, int m, int n)
Compute singular values when m >> n.
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_sv2uv(double *d, double **a, double **u, int m, double **v, int n)
Compute the singular value transformation when m >> n.
double b
#define x