GRASS Programmer's Manual  6.5.svn(2014)-r66266
sparse_matrix.c
Go to the documentation of this file.
1
2 /*****************************************************************************
3  *
4  * MODULE: Grass Gmath Library
5  * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6  * soerengebbert <at> gmx <dot> de
7  *
8  * PURPOSE: functions to manage linear equation systems
9  * part of the gmath library
10  *
11  * COPYRIGHT: (C) 2000 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 ** spmatrix, int rows)
122 {
123  int i;
124
125  if (spmatrix) {
126  for (i = 0; i < rows; i++)
127  G_math_free_spvector(spmatrix[i]);
128
129  G_free(spmatrix);
130  spmatrix = 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
207 G_math_spvector **G_math_A_to_Asp(double **A, int rows, double epsilon)
208 {
209  int i, j;
210
211  int nonull, count = 0;
212
213  G_math_spvector **Asp = NULL;
214
215  Asp = G_math_alloc_spmatrix(rows);
216
217 #pragma omp parallel for schedule (static) private(i, j, nonull, count)
218  for (i = 0; i < rows; i++) {
219  nonull = 0;
220  /*Count the number of non zero entries */
221  for (j = 0; j < rows; j++) {
222  if (A[i][j] > epsilon)
223  nonull++;
224  }
225  /*Allocate the sparse vector and insert values */
226  G_math_spvector *v = G_math_alloc_spvector(nonull);
227
228  count = 0;
229  for (j = 0; j < rows; j++) {
230  if (A[i][j] > epsilon) {
231  v->index[count] = j;
232  v->values[count] = A[i][j];
233  count++;
234  }
235  }
236  /*Add vector to sparse matrix */
238  }
239  return Asp;
240 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
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
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
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 **spmatrix, int rows)
Release the memory of the sparse matrix.
return NULL
Definition: dbfopen.c:1394
tuple cols
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse linear equation system 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.