GRASS Programmer's Manual  6.4.2(2012)
N_les.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines