GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
matrix.c
Go to the documentation of this file.
1 /*!
2  * \author
3  * Lubos Mitas (original program and various modifications)
4  *
5  * \author
6  * H. Mitasova,
7  * I. Kosinovsky, D. Gerdes,
8  * D. McCauley
9  * (GRASS4.1 version of the program and GRASS4.2 modifications)
10  *
11  * \author
12  * L. Mitas,
13  * H. Mitasova,
14  * I. Kosinovsky,
15  * D.Gerdes,
16  * D. McCauley
17  * (1993, 1995)
18  *
19  * \author modified by McCauley in August 1995
20  * \author modified by Mitasova in August 1995, Nov. 1996
21  *
22  * \copyright
23  * (C) 1993-1996 by Lubos Mitas and the GRASS Development Team
24  *
25  * \copyright
26  * This program is free software under the GNU General Public License (>=v2).
27  * Read the file COPYING that comes with GRASS for details.
28  */
29 
30 
31 #include <stdio.h>
32 #include <math.h>
33 #include <unistd.h>
34 #include <grass/gis.h>
35 #include <grass/interpf.h>
36 #include <grass/gmath.h>
37 
38 
39 int IL_matrix_create(struct interp_params *params,
40  struct triple *points, /* points for interpolation */
41  int n_points, /* number of points */
42  double **matrix, /* matrix */
43  int *indx)
44 {
45  static double *A = NULL;
46 
47  if (!A) {
48  if (!
49  (A =
50  G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
51  fprintf(stderr, "Cannot allocate memory for A\n");
52  return -1;
53  }
54  }
55  return IL_matrix_create_alloc(params, points, n_points, matrix, indx, A);
56 }
57 
58 /*!
59  * \brief Creates system of linear equations from interpolated points
60  *
61  * Creates system of linear equations represented by matrix using given
62  * points and interpolating function interp()
63  *
64  * \param params struct interp_params *
65  * \param points points for interpolation as struct triple
66  * \param n_points number of points
67  * \param[out] matrix the matrix
68  * \param indx
69  *
70  * \return -1 on failure, 1 on success
71  */
73  struct triple *points, /* points for interpolation */
74  int n_points, /* number of points */
75  double **matrix, /* matrix */
76  int *indx,
77  double *A /* temporary matrix unique for all threads */)
78 {
79  double xx, yy;
80  double rfsta2, r;
81  double d;
82  int n1, k1, k2, k, i1, l, m, i, j;
83  double fstar2 = params->fi * params->fi / 4.;
84  double RO, amaxa;
85  double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
86  double xxr, yyr;
87 
88  if (params->theta) {
89  teta = params->theta * (M_PI / 180); /* deg to rad */
90  rsin = sin(teta);
91  rcos = cos(teta);
92  }
93  if (params->scalex)
94  scale = params->scalex;
95 
96 
97  n1 = n_points + 1;
98 
99  /*
100  C GENERATION OF MATRIX
101  C FIRST COLUMN
102  */
103  A[1] = 0.;
104  for (k = 1; k <= n_points; k++) {
105  i1 = k + 1;
106  A[i1] = 1.;
107  }
108  /*
109  C OTHER COLUMNS
110  */
111  RO = -params->rsm;
112  /* fprintf (stderr, "sm[%d] = %f, ro=%f\n", 1, points[1].smooth, RO); */
113  for (k = 1; k <= n_points; k++) {
114  k1 = k * n1 + 1;
115  k2 = k + 1;
116  i1 = k1 + k;
117  if (params->rsm < 0.) { /*indicates variable smoothing */
118  A[i1] = -points[k - 1].sm; /* added by Mitasova nov. 96 */
119  /* G_debug(5, "sm[%d]=%f, a=%f", k, points[k-1].sm, A[i1]); */
120  }
121  else {
122  A[i1] = RO; /* constant smoothing */
123  }
124  /* if (i1 == 100) fprintf (stderr,i "A[%d] = %f\n", i1, A[i1]); */
125 
126  /* A[i1] = RO; */
127  for (l = k2; l <= n_points; l++) {
128  xx = points[k - 1].x - points[l - 1].x;
129  yy = points[k - 1].y - points[l - 1].y;
130 
131  if ((params->theta) && (params->scalex)) {
132  /* re run anisotropy */
133  xxr = xx * rcos + yy * rsin;
134  yyr = yy * rcos - xx * rsin;
135  xx = xxr;
136  yy = yyr;
137  r = scale * xx * xx + yy * yy;
138  rfsta2 = fstar2 * (scale * xx * xx + yy * yy);
139  }
140  else {
141  r = xx * xx + yy * yy;
142  rfsta2 = fstar2 * (xx * xx + yy * yy);
143  }
144 
145  if (rfsta2 == 0.) {
146  fprintf(stderr, "ident. points in segm.\n");
147  fprintf(stderr, "x[%d]=%f, x[%d]=%f, y[%d]=%f, y[%d]=%f\n",
148  k - 1, points[k - 1].x, l - 1, points[l - 1].x, k - 1,
149  points[k - 1].y, l - 1, points[l - 1].y);
150  return -1;
151  }
152  i1 = k1 + l;
153  A[i1] = params->interp(r, params->fi);
154  }
155  }
156 
157  /* C SYMMETRISATION */
158  amaxa = 1.;
159  for (k = 1; k <= n1; k++) {
160  k1 = (k - 1) * n1;
161  k2 = k + 1;
162  for (l = k2; l <= n1; l++) {
163  m = (l - 1) * n1 + k;
164  A[m] = A[k1 + l];
165  amaxa = amax1(A[m], amaxa);
166  }
167  }
168  m = 0;
169  for (i = 0; i <= n_points; i++) {
170  for (j = 0; j <= n_points; j++) {
171  m++;
172  matrix[i][j] = A[m];
173  }
174  }
175 
176  G_debug(3, "calling G_ludcmp() n=%d indx=%d", n_points, *indx);
177  if (G_ludcmp(matrix, n_points + 1, indx, &d) <= 0) {
178  /* find the inverse of the matrix */
179  fprintf(stderr, "G_ludcmp() failed! n=%d d=%.2f\n", n_points, d);
180  return -1;
181  }
182 
183  /* G_free_vector(A); */
184  return 1;
185 }
interp_fn * interp
Definition: interpf.h:101
double rsm
Definition: interpf.h:83
double scalex
Definition: interpf.h:90
double theta
Definition: interpf.h:89
int IL_matrix_create(struct interp_params *params, struct triple *points, int n_points, double **matrix, int *indx)
Definition: matrix.c:39
double amax1(double, double)
Definition: minmax.c:54
#define M_PI
Definition: gis.h:144
#define NULL
Definition: ccmath.h:32
#define x
double l
Definition: r_raster.c:39
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:41
double x
Definition: dataquad.h:42
double sm
Definition: dataquad.h:45
double fi
Definition: interpf.h:80
int G_ludcmp(double **, int, int *, double *)
LU decomposition.
Definition: lu.c:18
int IL_matrix_create_alloc(struct interp_params *params, struct triple *points, int n_points, double **matrix, int *indx, double *A)
Creates system of linear equations from interpolated points.
Definition: matrix.c:72
int G_debug(int, const char *,...) __attribute__((format(printf
double y
Definition: dataquad.h:43
double r
Definition: r_raster.c:39