Graphviz  2.41.20170921.2350
solve.c
Go to the documentation of this file.
1 /* $Id$ $Revision$ */
2 /* vim:set shiftwidth=4 ts=8: */
3 
4 /*************************************************************************
5  * Copyright (c) 2011 AT&T Intellectual Property
6  * All rights reserved. This program and the accompanying materials
7  * are made available under the terms of the Eclipse Public License v1.0
8  * which accompanies this distribution, and is available at
9  * http://www.eclipse.org/legal/epl-v10.html
10  *
11  * Contributors: See CVS logs. Details at http://www.graphviz.org/
12  *************************************************************************/
13 
14 
15 /* solves the system ab=c using gauss reduction */
16 #include <math.h>
17 #include <stdlib.h>
18 #include <stdio.h>
19 #include "render.h"
20 #define asub(i,j) a[(i)*n + (j)]
21 
22 
23 void solve(double *a, double *b, double *c, int n)
24 { /*a[n][n],b[n],c[n] */
25  double *asave, *csave;
26  double amax, dum, pivot;
27  register int i, ii, j;
28  register int k, m, mp;
29  register int istar, ip;
30  register int nm, nsq, t;
31 
32  istar = 0; /* quiet warnings */
33  nsq = n * n;
34  asave = N_GNEW(nsq, double);
35  csave = N_GNEW(n, double);
36 
37  for (i = 0; i < n; i++)
38  csave[i] = c[i];
39  for (i = 0; i < nsq; i++)
40  asave[i] = a[i];
41  /* eliminate ith unknown */
42  nm = n - 1;
43  for (i = 0; i < nm; i++) {
44  /* find largest pivot */
45  amax = 0.;
46  for (ii = i; ii < n; ii++) {
47  dum = fabs(asub(ii, i));
48  if (dum < amax)
49  continue;
50  istar = ii;
51  amax = dum;
52  }
53  /* return if pivot is too small */
54  if (amax < 1.e-10)
55  goto bad;
56  /* switch rows */
57  for (j = i; j < n; j++) {
58  t = istar * n + j;
59  dum = a[t];
60  a[t] = a[i * n + j];
61  a[i * n + j] = dum;
62  }
63  dum = c[istar];
64  c[istar] = c[i];
65  c[i] = dum;
66  /*pivot */
67  ip = i + 1;
68  for (ii = ip; ii < n; ii++) {
69  pivot = a[ii * n + i] / a[i * n + i];
70  c[ii] = c[ii] - pivot * c[i];
71  for (j = 0; j < n; j++)
72  a[ii * n + j] = a[ii * n + j] - pivot * a[i * n + j];
73  }
74  }
75  /* return if last pivot is too small */
76  if (fabs(a[n * n - 1]) < 1.e-10)
77  goto bad;
78  b[n - 1] = c[n - 1] / a[n * n - 1];
79  /* back substitute */
80  for (k = 0; k < nm; k++) {
81  m = n - k - 2;
82  b[m] = c[m];
83  mp = m + 1;
84  for (j = mp; j < n; j++)
85  b[m] = b[m] - a[m * n + j] * b[j];
86  b[m] = b[m] / a[m * n + m];
87  }
88  /* restore original a,c */
89  for (i = 0; i < n; i++)
90  c[i] = csave[i];
91  for (i = 0; i < nsq; i++)
92  a[i] = asave[i];
93  free(asave);
94  free(csave);
95  return;
96  bad:
97  printf("ill-conditioned\n");
98  free(asave);
99  free(csave);
100  return;
101 }
void solve(double *, double *, double *, int)
Definition: solve.c:23
#define asub(i, j)
Definition: solve.c:20
#define N_GNEW(n, t)
Definition: agxbuf.c:20