GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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  * */
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  * */
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  * */
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  * */
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  * */
121 void 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  * */
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  * */
179 double **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  * */
221 double **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  * */
253 G_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 */
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  * */
304 G_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 
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  * */
363 void 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 }
The row vector of the sparse matrix.
Definition: gmath.h:54
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_alloc_spvector(int cols)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:76
G_math_spvector ** G_math_A_to_Asp(double **A, int rows, double epsilon)
Convert a quadratic matrix into a sparse matrix.
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
Definition: dalloc.c:60
#define NULL
Definition: ccmath.h:32
#define x
#define G_calloc(m, n)
Definition: defs/gis.h:113
unsigned int * index
Definition: gmath.h:58
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
double * values
Definition: gmath.h:56
unsigned int cols
Definition: gmath.h:57
void G_math_free_spvector(G_math_spvector *spvector)
Release the memory of the sparse vector.
Definition: sparse_matrix.c:98
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_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_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
int G_debug(int, const char *,...) __attribute__((format(printf
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.