GRASS Programmer's Manual  6.4.2(2012)
N_les_pivot.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:      linear equation system pivoting strategy
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 <math.h>
00020 #include <unistd.h>
00021 #include <stdio.h>
00022 #include <string.h>
00023 #include "grass/N_pde.h"
00024 #include "solvers_local_proto.h"
00025 
00026 
00027 #define TINY 1.0e-20
00028 
00029 
00030 
00047 int N_les_pivot_create(N_les * les)
00048 {
00049     int num = 0;                /*number of changed rows */
00050     int i, j, k;
00051     double max;
00052     int number = 0;
00053     double tmpval = 0.0, s = 0.0;
00054     double *link = NULL;
00055 
00056     G_debug(2, "N_les_pivot_create: swap rows if needed");
00057     for (i = 0; i < les->rows; i++) {
00058         max = fabs(les->A[i][i]);
00059         number = i;
00060         for (j = i; j < les->rows; j++) {
00061             s = 0.0;
00062             for (k = i; k < les->rows; k++) {
00063                 s += fabs(les->A[j][i]);
00064             }
00065             /*search for the pivot element */
00066             if (max < fabs(les->A[j][i]) / s) {
00067                 max = fabs(les->A[j][i]);
00068                 number = j;
00069             }
00070         }
00071         if (max == 0) {
00072             max = TINY;
00073             G_warning("Matrix is singular");
00074         }
00075         /*if an pivot element was found, swap the les entries */
00076         if (number != i) {
00077 
00078             G_debug(4, "swap row %i with row %i", i, number);
00079 
00080             tmpval = les->b[number];
00081             les->b[number] = les->b[i];
00082             les->b[i] = tmpval;
00083 
00084             link = les->A[number];
00085             les->A[number] = les->A[i];
00086             les->A[i] = link;
00087             num++;
00088         }
00089     }
00090 
00091     return num;
00092 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines