GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /***************************************************************************** 00003 * 00004 * MODULE: Grass PDE Numerical Library 00005 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006 00006 * soerengebbert <at> gmx <dot> de 00007 * 00008 * PURPOSE: functions to manage linear equation systems 00009 * part of the gpde library 00010 * 00011 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 * 00017 *****************************************************************************/ 00018 00019 #include "grass/N_pde.h" 00020 #include <stdlib.h> 00021 00029 N_spvector *N_alloc_spvector(int cols) 00030 { 00031 N_spvector *spvector; 00032 00033 G_debug(4, "Allocate memory for a sparse vector with %i cols\n", cols); 00034 00035 spvector = (N_spvector *) G_calloc(1, sizeof(N_spvector)); 00036 00037 spvector->cols = cols; 00038 spvector->index = (int *)G_calloc(cols, sizeof(int)); 00039 spvector->values = (double *)G_calloc(cols, sizeof(double)); 00040 00041 return spvector; 00042 } 00043 00055 N_les *N_alloc_nquad_les(int cols, int rows, int type) 00056 { 00057 return N_alloc_les_param(cols, rows, type, 2); 00058 } 00059 00071 N_les *N_alloc_nquad_les_Ax(int cols, int rows, int type) 00072 { 00073 return N_alloc_les_param(cols, rows, type, 1); 00074 } 00075 00087 N_les *N_alloc_nquad_les_A(int cols, int rows, int type) 00088 { 00089 return N_alloc_les_param(cols, rows, type, 0); 00090 } 00091 00103 N_les *N_alloc_nquad_les_Ax_b(int cols, int rows, int type) 00104 { 00105 return N_alloc_les_param(cols, rows, type, 2); 00106 } 00107 00108 00109 00120 N_les *N_alloc_les(int rows, int type) 00121 { 00122 return N_alloc_les_param(rows, rows, type, 2); 00123 } 00124 00135 N_les *N_alloc_les_Ax(int rows, int type) 00136 { 00137 return N_alloc_les_param(rows, rows, type, 1); 00138 } 00139 00150 N_les *N_alloc_les_A(int rows, int type) 00151 { 00152 return N_alloc_les_param(rows, rows, type, 0); 00153 } 00154 00165 N_les *N_alloc_les_Ax_b(int rows, int type) 00166 { 00167 return N_alloc_les_param(rows, rows, type, 2); 00168 } 00169 00170 00198 N_les *N_alloc_les_param(int cols, int rows, int type, int parts) 00199 { 00200 N_les *les; 00201 int i; 00202 00203 if (type == N_SPARSE_LES) 00204 G_debug(2, 00205 "Allocate memory for a sparse linear equation system with %i rows\n", 00206 rows); 00207 else 00208 G_debug(2, 00209 "Allocate memory for a regular linear equation system with %i rows\n", 00210 rows); 00211 00212 les = (N_les *) G_calloc(1, sizeof(N_les)); 00213 00214 if (parts > 0) { 00215 les->x = (double *)G_calloc(cols, sizeof(double)); 00216 for (i = 0; i < cols; i++) 00217 les->x[i] = 0.0; 00218 } 00219 00220 00221 if (parts > 1) { 00222 les->b = (double *)G_calloc(cols, sizeof(double)); 00223 for (i = 0; i < cols; i++) 00224 les->b[i] = 0.0; 00225 } 00226 00227 les->A = NULL; 00228 les->Asp = NULL; 00229 les->rows = rows; 00230 les->cols = cols; 00231 if (rows == cols) 00232 les->quad = 1; 00233 else 00234 les->quad = 0; 00235 00236 if (type == N_SPARSE_LES) { 00237 les->Asp = (N_spvector **) G_calloc(rows, sizeof(N_spvector *)); 00238 les->type = N_SPARSE_LES; 00239 } 00240 else { 00241 les->A = (double **)G_calloc(rows, sizeof(double *)); 00242 for (i = 0; i < rows; i++) { 00243 les->A[i] = (double *)G_calloc(cols, sizeof(double)); 00244 } 00245 les->type = N_NORMAL_LES; 00246 } 00247 00248 return les; 00249 } 00250 00251 00263 int N_add_spvector_to_les(N_les * les, N_spvector * spvector, int row) 00264 { 00265 00266 00267 if (les != NULL) { 00268 if (les->type != N_SPARSE_LES) 00269 return -1; 00270 00271 if (les->rows > row) { 00272 G_debug(5, 00273 "Add sparse vector %p to the sparse linear equation system at row %i\n", 00274 spvector, row); 00275 les->Asp[row] = spvector; 00276 } 00277 else 00278 return -1; 00279 00280 } 00281 else { 00282 return -1; 00283 } 00284 00285 00286 return 1; 00287 } 00288 00312 void N_print_les(N_les * les) 00313 { 00314 int i, j, k, out; 00315 00316 00317 if (les->type == N_SPARSE_LES) { 00318 for (i = 0; i < les->rows; i++) { 00319 for (j = 0; j < les->cols; j++) { 00320 out = 0; 00321 for (k = 0; k < les->Asp[i]->cols; k++) { 00322 if (les->Asp[i]->index[k] == j) { 00323 fprintf(stdout, "%4.5f ", les->Asp[i]->values[k]); 00324 out = 1; 00325 } 00326 } 00327 if (!out) 00328 fprintf(stdout, "%4.5f ", 0.0); 00329 } 00330 if (les->x) 00331 fprintf(stdout, " * %4.5f", les->x[i]); 00332 if (les->b) 00333 fprintf(stdout, " = %4.5f ", les->b[i]); 00334 00335 fprintf(stdout, "\n"); 00336 } 00337 } 00338 else { 00339 00340 for (i = 0; i < les->rows; i++) { 00341 for (j = 0; j < les->cols; j++) { 00342 fprintf(stdout, "%4.5f ", les->A[i][j]); 00343 } 00344 if (les->x) 00345 fprintf(stdout, " * %4.5f", les->x[i]); 00346 if (les->b) 00347 fprintf(stdout, " = %4.5f ", les->b[i]); 00348 00349 fprintf(stdout, "\n"); 00350 } 00351 00352 } 00353 return; 00354 } 00355 00363 void N_free_spvector(N_spvector * spvector) 00364 { 00365 if (spvector) { 00366 if (spvector->values) 00367 G_free(spvector->values); 00368 if (spvector->index) 00369 G_free(spvector->index); 00370 G_free(spvector); 00371 00372 spvector = NULL; 00373 } 00374 00375 return; 00376 } 00377 00378 00387 void N_free_les(N_les * les) 00388 { 00389 int i; 00390 00391 if (les->type == N_SPARSE_LES) 00392 G_debug(2, "Releasing memory of a sparse linear equation system\n"); 00393 else 00394 G_debug(2, "Releasing memory of a regular linear equation system\n"); 00395 00396 if (les) { 00397 00398 if (les->x) 00399 G_free(les->x); 00400 if (les->b) 00401 G_free(les->b); 00402 00403 if (les->type == N_SPARSE_LES) { 00404 00405 if (les->Asp) { 00406 for (i = 0; i < les->rows; i++) 00407 if (les->Asp[i]) 00408 N_free_spvector(les->Asp[i]); 00409 00410 G_free(les->Asp); 00411 } 00412 } 00413 else { 00414 00415 if (les->A) { 00416 for (i = 0; i < les->rows; i++) 00417 if (les->A[i]) 00418 G_free(les->A[i]); 00419 00420 G_free(les->A); 00421 } 00422 } 00423 00424 free(les); 00425 } 00426 00427 return; 00428 }