GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
sparse_matrix.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: linear equation system solvers
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#include <stdlib.h>
20#include <math.h>
21#include <grass/gmath.h>
22#include <grass/gis.h>
23
24/*!
25 * \brief Adds a sparse vector to a sparse matrix at position row
26 *
27 * Return 1 for success and -1 for failure
28 *
29 * \param Asp G_math_spvector **
30 * \param spvector G_math_spvector *
31 * \param row int
32 * \return int 1 success, -1 failure
33 *
34 * */
35int G_math_add_spvector(G_math_spvector ** Asp, G_math_spvector * spvector,
36 int row)
37{
38 if (Asp != NULL) {
39 G_debug(5,
40 "Add sparse vector %p to the sparse linear equation system at row %i\n",
41 spvector, row);
42 Asp[row] = spvector;
43 }
44 else {
45 return -1;
46 }
47
48 return 1;
49}
50
51/*!
52 * \brief Allocate memory for a sparse matrix
53 *
54 * \param rows int
55 * \return G_math_spvector **
56 *
57 * */
58G_math_spvector **G_math_alloc_spmatrix(int rows)
59{
60 G_math_spvector **spmatrix;
61
62 G_debug(4, "Allocate memory for a sparse matrix with %i rows\n", rows);
63
64 spmatrix = (G_math_spvector **) G_calloc(rows, sizeof(G_math_spvector *));
65
66 return spmatrix;
67}
68
69/*!
70 * \brief Allocate memory for a sparse vector
71 *
72 * \param cols int
73 * \return G_math_spvector *
74 *
75 * */
76G_math_spvector *G_math_alloc_spvector(int cols)
77{
78 G_math_spvector *spvector;
79
80 G_debug(4, "Allocate memory for a sparse vector with %i cols\n", cols);
81
82 spvector = (G_math_spvector *) G_calloc(1, sizeof(G_math_spvector));
83
84 spvector->cols = cols;
85 spvector->index = (unsigned int *)G_calloc(cols, sizeof(unsigned int));
86 spvector->values = (double *)G_calloc(cols, sizeof(double));
87
88 return spvector;
89}
90
91/*!
92 * \brief Release the memory of the sparse vector
93 *
94 * \param spvector G_math_spvector *
95 * \return void
96 *
97 * */
98void G_math_free_spvector(G_math_spvector * spvector)
99{
100 if (spvector) {
101 if (spvector->values)
102 G_free(spvector->values);
103 if (spvector->index)
104 G_free(spvector->index);
105 G_free(spvector);
106
107 spvector = NULL;
108 }
109
110 return;
111}
112
113/*!
114 * \brief Release the memory of the sparse matrix
115 *
116 * \param Asp G_math_spvector **
117 * \param rows int
118 * \return void
119 *
120 * */
121void G_math_free_spmatrix(G_math_spvector ** Asp, int rows)
122{
123 int i;
124
125 if (Asp) {
126 for (i = 0; i < rows; i++)
127 G_math_free_spvector(Asp[i]);
128
129 G_free(Asp);
130 Asp = NULL;
131 }
132
133 return;
134}
135
136/*!
137 *
138 * \brief print the sparse matrix Asp to stdout
139 *
140 *
141 * \param Asp (G_math_spvector **)
142 * \param rows (int)
143 * \return void
144 *
145 * */
146void G_math_print_spmatrix(G_math_spvector ** Asp, int rows)
147{
148 int i, j, k, out;
149
150 for (i = 0; i < rows; i++) {
151 for (j = 0; j < rows; j++) {
152 out = 0;
153 for (k = 0; k < Asp[i]->cols; k++) {
154 if (Asp[i]->index[k] == j) {
155 fprintf(stdout, "%4.5f ", Asp[i]->values[k]);
156 out = 1;
157 }
158 }
159 if (!out)
160 fprintf(stdout, "%4.5f ", 0.0);
161 }
162 fprintf(stdout, "\n");
163 }
164
165 return;
166}
167
168
169/*!
170 * \brief Convert a sparse matrix into a quadratic matrix
171 *
172 * This function is multi-threaded with OpenMP. It creates its own parallel OpenMP region.
173 *
174 * \param Asp (G_math_spvector **)
175 * \param rows (int)
176 * \return (double **)
177 *
178 * */
179double **G_math_Asp_to_A(G_math_spvector ** Asp, int rows)
180{
181 int i, j;
182
183 double **A = NULL;
184
185 A = G_alloc_matrix(rows, rows);
186
187#pragma omp parallel for schedule (static) private(i, j)
188 for (i = 0; i < rows; i++) {
189 for (j = 0; j < Asp[i]->cols; j++) {
190 A[i][Asp[i]->index[j]] = Asp[i]->values[j];
191 }
192 }
193 return A;
194}
195
196/*!
197 * \brief Convert a symmetric sparse matrix into a symmetric band matrix
198 *
199 \verbatim
200 Symmetric matrix with bandwidth of 3
201
202 5 2 1 0
203 2 5 2 1
204 1 2 5 2
205 0 1 2 5
206
207 will be converted into the band matrix
208
209 5 2 1
210 5 2 1
211 5 2 0
212 5 0 0
213
214 \endverbatim
215 * \param Asp (G_math_spvector **)
216 * \param rows (int)
217 * \param bandwidth (int)
218 * \return (double **) the resulting ymmetric band matrix [rows][bandwidth]
219 *
220 * */
221double **G_math_Asp_to_sband_matrix(G_math_spvector ** Asp, int rows, int bandwidth)
222{
223 int i, j;
224
225 double **A = NULL;
226
227 A = G_alloc_matrix(rows, bandwidth);
228
229 for (i = 0; i < rows; i++) {
230 for (j = 0; j < Asp[i]->cols; j++) {
231 if(Asp[i]->index[j] == i) {
232 A[i][0] = Asp[i]->values[j];
233 } else if (Asp[i]->index[j] > i) {
234 A[i][Asp[i]->index[j] - i] = Asp[i]->values[j];
235 }
236 }
237 }
238 return A;
239}
240
241
242/*!
243 * \brief Convert a quadratic matrix into a sparse matrix
244 *
245 * This function is multi-threaded with OpenMP. It creates its own parallel OpenMP region.
246 *
247 * \param A (double **)
248 * \param rows (int)
249 * \param epsilon (double) -- non-zero values are greater then epsilon
250 * \return (G_math_spvector **)
251 *
252 * */
253G_math_spvector **G_math_A_to_Asp(double **A, int rows, double epsilon)
254{
255 int i, j;
256
257 int nonull, count = 0;
258
259 G_math_spvector **Asp = NULL;
260
261 Asp = G_math_alloc_spmatrix(rows);
262
263#pragma omp parallel for schedule (static) private(i, j, nonull, count)
264 for (i = 0; i < rows; i++) {
265 nonull = 0;
266 /*Count the number of non zero entries */
267 for (j = 0; j < rows; j++) {
268 if (A[i][j] > epsilon)
269 nonull++;
270 }
271 /*Allocate the sparse vector and insert values */
272 G_math_spvector *v = G_math_alloc_spvector(nonull);
273
274 count = 0;
275 for (j = 0; j < rows; j++) {
276 if (A[i][j] > epsilon) {
277 v->index[count] = j;
278 v->values[count] = A[i][j];
279 count++;
280 }
281 }
282 /*Add vector to sparse matrix */
283 G_math_add_spvector(Asp, v, i);
284 }
285 return Asp;
286}
287
288
289
290/*!
291 * \brief Convert a symmetric band matrix into a sparse matrix
292 *
293 * WARNING:
294 * This function is experimental, do not use.
295 * Only the upper triangle matrix of the band strcuture is copied.
296 *
297 * \param A (double **) the symmetric band matrix
298 * \param rows (int)
299 * \param bandwidth (int)
300 * \param epsilon (double) -- non-zero values are greater then epsilon
301 * \return (G_math_spvector **)
302 *
303 * */
304G_math_spvector **G_math_sband_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)
305{
306 int i, j;
307
308 int nonull, count = 0;
309
310 G_math_spvector **Asp = NULL;
311
312 Asp = G_math_alloc_spmatrix(rows);
313
314 for (i = 0; i < rows; i++) {
315 nonull = 0;
316 /*Count the number of non zero entries */
317 for (j = 0; j < bandwidth; j++) {
318 if (A[i][j] > epsilon)
319 nonull++;
320 }
321
322 /*Allocate the sparse vector and insert values */
323
324 G_math_spvector *v = G_math_alloc_spvector(nonull);
325
326 count = 0;
327 if (A[i][0] > epsilon) {
328 v->index[count] = i;
329 v->values[count] = A[i][0];
330 count++;
331 }
332
333 for (j = 1; j < bandwidth; j++) {
334 if (A[i][j] > epsilon && i + j < rows) {
335 v->index[count] = i + j;
336 v->values[count] = A[i][j];
337 count++;
338 }
339 }
340 /*Add vector to sparse matrix */
341 G_math_add_spvector(Asp, v, i);
342 }
343 return Asp;
344}
345
346
347/*!
348 * \brief Compute the matrix - vector product
349 * of sparse matrix **Asp and vector x.
350 *
351 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
352 *
353 * y = A * x
354 *
355 *
356 * \param Asp (G_math_spvector **)
357 * \param x (double) *)
358 * \param y (double * )
359 * \param rows (int)
360 * \return (void)
361 *
362 * */
363void G_math_Ax_sparse(G_math_spvector ** Asp, double *x, double *y, int rows)
364{
365 int i, j;
366
367 double tmp;
368
369#pragma omp for schedule (static) private(i, j, tmp)
370 for (i = 0; i < rows; i++) {
371 tmp = 0;
372 for (j = 0; j < Asp[i]->cols; j++) {
373 tmp += Asp[i]->values[j] * x[Asp[i]->index[j]];
374 }
375 y[i] = tmp;
376 }
377 return;
378}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#define NULL
Definition: ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
int count
void G_math_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
double ** G_math_Asp_to_sband_matrix(G_math_spvector **Asp, int rows, int bandwidth)
Convert a symmetric sparse matrix into a symmetric band matrix.
void G_math_print_spmatrix(G_math_spvector **Asp, int rows)
print the sparse matrix Asp to stdout
double ** G_math_Asp_to_A(G_math_spvector **Asp, int rows)
Convert a sparse matrix into a quadratic matrix.
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:76
void G_math_free_spvector(G_math_spvector *spvector)
Release the memory of the sparse vector.
Definition: sparse_matrix.c:98
void G_math_free_spmatrix(G_math_spvector **Asp, int rows)
Release the memory of the sparse matrix.
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse matrix at position row.
Definition: sparse_matrix.c:35
G_math_spvector ** G_math_sband_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)
Convert a symmetric band matrix into a sparse matrix.
G_math_spvector ** G_math_A_to_Asp(double **A, int rows, double epsilon)
Convert a quadratic matrix into a sparse matrix.
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
#define x