GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
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
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}
#define NULL
Definition: ccmath.h:32
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double l
double r
double amax1(double, double)
Definition: minmax.c:54
int G_ludcmp(double **a, int n, int *indx, double *d)
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 IL_matrix_create(struct interp_params *params, struct triple *points, int n_points, double **matrix, int *indx)
Definition: matrix.c:39
interp_fn * interp
Definition: interpf.h:101
double fi
Definition: interpf.h:80
double theta
Definition: interpf.h:89
double rsm
Definition: interpf.h:83
double scalex
Definition: interpf.h:90
double sm
Definition: dataquad.h:45
double x
Definition: dataquad.h:42
double y
Definition: dataquad.h:43
#define x