GRASS GIS 7 Programmer's Manual  7.0.5(2016)-r00000
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 
184 int G_math_solv(double **a,double *b,int n)
185 {
186  return solv(a[0],b, n);
187 }
188 
189 
198  int G_math_solvps(double **a,double *b,int n)
199 {
200  return solvps(a[0], b,n);
201 }
202 
203 
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 
242 int G_math_minv(double **a,int n)
243 {
244  return minv(a[0], n);
245 }
246 
247 
256 int G_math_psinv(double **a,int n)
257 {
258  return psinv( a[0], n);
259 }
260 
261 
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 */
289 void G_math_eigval(double **a,double *ev,int n)
290 {
291  eigval(a[0], ev, n);
292  return;
293 }
294 
295 
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 
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 
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 
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 
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.
double b
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
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