1
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical 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 gpde 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 <grass/N_pde.h>
21 #include <grass/gmath.h>
22
23
24 /*!
25 * \brief Allocate memory for a (not) quadratic linear equation system which includes the Matrix A, vector x and vector b
26 *
27 * This function calls #N_alloc_les_param
28 *
29 * \param cols int
30 * \param rows int
31 * \param type int
32 * \return N_les *
33 *
34 * */
N_alloc_nquad_les(int cols,int rows,int type)35 N_les *N_alloc_nquad_les(int cols, int rows, int type)
36 {
37 return N_alloc_les_param(cols, rows, type, 2);
38 }
39
40 /*!
41 * \brief Allocate memory for a (not) quadratic linear equation system which includes the Matrix A and vector x
42 *
43 * This function calls #N_alloc_les_param
44 *
45 * \param cols int
46 * \param rows int
47 * \param type int
48 * \return N_les *
49 *
50 * */
N_alloc_nquad_les_Ax(int cols,int rows,int type)51 N_les *N_alloc_nquad_les_Ax(int cols, int rows, int type)
52 {
53 return N_alloc_les_param(cols, rows, type, 1);
54 }
55
56 /*!
57 * \brief Allocate memory for a (not) quadratic linear equation system which includes the Matrix A
58 *
59 * This function calls #N_alloc_les_param
60 *
61 * \param cols int
62 * \param rows int
63 * \param type int
64 * \return N_les *
65 *
66 * */
N_alloc_nquad_les_A(int cols,int rows,int type)67 N_les *N_alloc_nquad_les_A(int cols, int rows, int type)
68 {
69 return N_alloc_les_param(cols, rows, type, 0);
70 }
71
72 /*!
73 * \brief Allocate memory for a (not) quadratic linear equation system which includes the Matrix A, vector x and vector b
74 *
75 * This function calls #N_alloc_les_param
76 *
77 * \param cols int
78 * \param rows int
79 * \param type int
80 * \return N_les *
81 *
82 * */
N_alloc_nquad_les_Ax_b(int cols,int rows,int type)83 N_les *N_alloc_nquad_les_Ax_b(int cols, int rows, int type)
84 {
85 return N_alloc_les_param(cols, rows, type, 2);
86 }
87
88
89
90 /*!
91 * \brief Allocate memory for a quadratic linear equation system which includes the Matrix A, vector x and vector b
92 *
93 * This function calls #N_alloc_les_param
94 *
95 * \param rows int
96 * \param type int
97 * \return N_les *
98 *
99 * */
N_alloc_les(int rows,int type)100 N_les *N_alloc_les(int rows, int type)
101 {
102 return N_alloc_les_param(rows, rows, type, 2);
103 }
104
105 /*!
106 * \brief Allocate memory for a quadratic linear equation system which includes the Matrix A and vector x
107 *
108 * This function calls #N_alloc_les_param
109 *
110 * \param rows int
111 * \param type int
112 * \return N_les *
113 *
114 * */
N_alloc_les_Ax(int rows,int type)115 N_les *N_alloc_les_Ax(int rows, int type)
116 {
117 return N_alloc_les_param(rows, rows, type, 1);
118 }
119
120 /*!
121 * \brief Allocate memory for a quadratic linear equation system which includes the Matrix A
122 *
123 * This function calls #N_alloc_les_param
124 *
125 * \param rows int
126 * \param type int
127 * \return N_les *
128 *
129 * */
N_alloc_les_A(int rows,int type)130 N_les *N_alloc_les_A(int rows, int type)
131 {
132 return N_alloc_les_param(rows, rows, type, 0);
133 }
134
135 /*!
136 * \brief Allocate memory for a quadratic linear equation system which includes the Matrix A, vector x and vector b
137 *
138 * This function calls #N_alloc_les_param
139 *
140 * \param rows int
141 * \param type int
142 * \return N_les *
143 *
144 * */
N_alloc_les_Ax_b(int rows,int type)145 N_les *N_alloc_les_Ax_b(int rows, int type)
146 {
147 return N_alloc_les_param(rows, rows, type, 2);
148 }
149
150
151 /*!
152 * \brief Allocate memory for a quadratic or not quadratic linear equation system
153 *
154 * The type of the linear equation system must be N_NORMAL_LES for
155 * a regular quadratic matrix or N_SPARSE_LES for a sparse matrix
156 *
157 * <p>
158 * In case of N_NORMAL_LES
159 *
160 * A quadratic matrix of size rows*rows*sizeof(double) will allocated
161 *
162 * <p>
163 * In case of N_SPARSE_LES
164 *
165 * a vector of size row will be allocated, ready to hold additional allocated sparse vectors.
166 * each sparse vector may have a different size.
167 *
168 * Parameter parts defines which parts of the les should be allocated.
169 * The number of columns and rows defines if the matrix is quadratic.
170 *
171 * \param cols int
172 * \param rows int
173 * \param type int
174 * \param parts int -- 2 = A, x and b; 1 = A and x; 0 = A allocated
175 * \return N_les *
176 *
177 * */
N_alloc_les_param(int cols,int rows,int type,int parts)178 N_les *N_alloc_les_param(int cols, int rows, int type, int parts)
179 {
180 N_les *les;
181
182 int i;
183
184 if (type == N_SPARSE_LES)
185 G_debug(2,
186 "Allocate memory for a sparse linear equation system with %i rows\n",
187 rows);
188 else
189 G_debug(2,
190 "Allocate memory for a regular linear equation system with %i rows\n",
191 rows);
192
193 les = (N_les *) G_calloc(1, sizeof(N_les));
194
195 if (parts > 0) {
196 les->x = (double *)G_calloc(cols, sizeof(double));
197 for (i = 0; i < cols; i++)
198 les->x[i] = 0.0;
199 }
200
201
202 if (parts > 1) {
203 les->b = (double *)G_calloc(cols, sizeof(double));
204 for (i = 0; i < cols; i++)
205 les->b[i] = 0.0;
206 }
207
208 les->A = NULL;
209 les->Asp = NULL;
210 les->rows = rows;
211 les->cols = cols;
212 if (rows == cols)
213 les->quad = 1;
214 else
215 les->quad = 0;
216
217 if (type == N_SPARSE_LES) {
218 les->Asp = G_math_alloc_spmatrix(rows);
219 les->type = N_SPARSE_LES;
220 }
221 else {
222 les->A = G_alloc_matrix(rows, cols);
223 les->type = N_NORMAL_LES;
224 }
225
226 return les;
227 }
228
229 /*!
230 *
231 * \brief prints the linear equation system to stdout
232 *
233 * <p>
234 * Format:
235 * A*x = b
236 *
237 * <p>
238 * Example
239 \verbatim
240
241 2 1 1 1 * 2 = 0.1
242 1 2 0 0 * 3 = 0.2
243 1 0 2 0 * 3 = 0.2
244 1 0 0 2 * 2 = 0.1
245
246 \endverbatim
247 *
248 * \param les N_les *
249 * \return void
250 *
251 * */
N_print_les(N_les * les)252 void N_print_les(N_les * les)
253 {
254 int i, j, k, out;
255
256
257 if (les->type == N_SPARSE_LES) {
258 for (i = 0; i < les->rows; i++) {
259 for (j = 0; j < les->cols; j++) {
260 out = 0;
261 for (k = 0; k < les->Asp[i]->cols; k++) {
262 if (les->Asp[i]->index[k] == j) {
263 fprintf(stdout, "%4.5f ", les->Asp[i]->values[k]);
264 out = 1;
265 }
266 }
267 if (!out)
268 fprintf(stdout, "%4.5f ", 0.0);
269 }
270 if (les->x)
271 fprintf(stdout, " * %4.5f", les->x[i]);
272 if (les->b)
273 fprintf(stdout, " = %4.5f ", les->b[i]);
274
275 fprintf(stdout, "\n");
276 }
277 }
278 else {
279
280 for (i = 0; i < les->rows; i++) {
281 for (j = 0; j < les->cols; j++) {
282 fprintf(stdout, "%4.5f ", les->A[i][j]);
283 }
284 if (les->x)
285 fprintf(stdout, " * %4.5f", les->x[i]);
286 if (les->b)
287 fprintf(stdout, " = %4.5f ", les->b[i]);
288
289 fprintf(stdout, "\n");
290 }
291
292 }
293 return;
294 }
295
296 /*!
297 * \brief Release the memory of the linear equation system
298 *
299 * \param les N_les *
300 * \return void
301 *
302 * */
303
N_free_les(N_les * les)304 void N_free_les(N_les * les)
305 {
306 if (les->type == N_SPARSE_LES)
307 G_debug(2, "Releasing memory of a sparse linear equation system\n");
308 else
309 G_debug(2, "Releasing memory of a regular linear equation system\n");
310
311 if (les) {
312
313 if (les->x)
314 G_free(les->x);
315 if (les->b)
316 G_free(les->b);
317
318 if (les->type == N_SPARSE_LES) {
319
320 if (les->Asp) {
321 G_math_free_spmatrix(les->Asp, les->rows);
322 }
323 }
324 else {
325
326 if (les->A) {
327 G_free_matrix(les->A);
328 }
329 }
330
331 free(les);
332 }
333
334 return;
335 }
336