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