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