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
23void 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}
102