GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-55b2a2bccc
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  * */
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  * */
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  * */
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  * */
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  * */
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  * */
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  * */
181 double **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  * */
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  * */
260 G_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 */
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  * */
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 
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  * */
369 void 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 }
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
int G_debug(int, const char *,...) __attribute__((format(printf
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
Definition: dalloc.c:57
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.
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_spvector(int cols)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:77
double ** G_math_Asp_to_A(G_math_spvector **Asp, int rows)
Convert a sparse matrix into a quadratic matrix.
void G_math_print_spmatrix(G_math_spvector **Asp, int rows)
print the sparse matrix Asp to stdout
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.
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:59
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.
void G_math_free_spvector(G_math_spvector *spvector)
Release the memory of the sparse vector.
Definition: sparse_matrix.c:99
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
The row vector of the sparse matrix.
Definition: gmath.h:54
double * values
Definition: gmath.h:55
unsigned int cols
Definition: gmath.h:56
unsigned int * index
Definition: gmath.h:57
#define x