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: 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 }