GRASS Programmer's Manual  6.4.4(2014)-r
N_les_pivot.c
Go to the documentation of this file.
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: linear equation system pivoting strategy
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 <math.h>
20 #include <unistd.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include "grass/N_pde.h"
24 #include "solvers_local_proto.h"
25 
26 
27 #define TINY 1.0e-20
28 
29 
30 
48 {
49  int num = 0; /*number of changed rows */
50  int i, j, k;
51  double max;
52  int number = 0;
53  double tmpval = 0.0, s = 0.0;
54  double *link = NULL;
55 
56  G_debug(2, "N_les_pivot_create: swap rows if needed");
57  for (i = 0; i < les->rows; i++) {
58  max = fabs(les->A[i][i]);
59  number = i;
60  for (j = i; j < les->rows; j++) {
61  s = 0.0;
62  for (k = i; k < les->rows; k++) {
63  s += fabs(les->A[j][i]);
64  }
65  /*search for the pivot element */
66  if (max < fabs(les->A[j][i]) / s) {
67  max = fabs(les->A[j][i]);
68  number = j;
69  }
70  }
71  if (max == 0) {
72  max = TINY;
73  G_warning("Matrix is singular");
74  }
75  /*if an pivot element was found, swap the les entries */
76  if (number != i) {
77 
78  G_debug(4, "swap row %i with row %i", i, number);
79 
80  tmpval = les->b[number];
81  les->b[number] = les->b[i];
82  les->b[i] = tmpval;
83 
84  link = les->A[number];
85  les->A[number] = les->A[i];
86  les->A[i] = link;
87  num++;
88  }
89  }
90 
91  return num;
92 }
#define NULL
Definition: strings.c:26
int N_les_pivot_create(N_les *les)
Optimize the structure of the linear equation system with a common pivoting strategy.
Definition: N_les_pivot.c:47
#define max(x, y)
Definition: draw2.c:69
int G_warning(const char *msg,...)
Print a warning message to stderr.
double ** A
Definition: N_pde.h:100
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
double * b
Definition: N_pde.h:99
int rows
Definition: N_pde.h:102
#define TINY
Definition: N_les_pivot.c:27
The linear equation system (les) structure.
Definition: N_pde.h:96