← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gpde/n_les.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gpde/n_les.c
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 * */
 
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
 
 
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
}