GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
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 
35 int 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 
58 G_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 
76 G_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 
98 void 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 
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 
146 void 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 
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 
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 
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 */
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 
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 
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 
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 }
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
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
G_math_spvector ** G_math_A_to_Asp(double **A, int rows, double epsilon)
Convert a quadratic matrix into a sparse matrix.
int count
#define NULL
Definition: ccmath.h:32
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
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
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.