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 | * This code was (mostly) written by Ken Turkowski, who said: |
16 | * |
17 | * Oh, that. I wrote it in college the first time. It's open source - I think I |
18 | * posted it after seeing so many people solve equations by inverting matrices |
19 | * by computing minors naïvely. |
20 | * -Ken |
21 | * |
22 | * The views represented here are mine and are not necessarily shared by |
23 | * my employer. |
24 | Ken Turkowski turk@apple.com |
25 | Immersive Media Technologist http://www.worldserver.com/turk/ |
26 | Apple Computer, Inc. |
27 | 1 Infinite Loop, MS 302-3VR |
28 | Cupertino, CA 95014 |
29 | */ |
30 | |
31 | /* Matinv() inverts the matrix A using LU decomposition. Arguments: |
32 | * A - the (n x n) matrix to be inverted |
33 | * Ainv - the (n x n) inverted matrix |
34 | * n - the order of the matrices A and Ainv |
35 | */ |
36 | |
37 | #include <stdlib.h> |
38 | #include "render.h" |
39 | extern int lu_decompose(double **a, int n); |
40 | extern void lu_solve(double *x, double *b, int n); |
41 | |
42 | int matinv(double **A, double **Ainv, int n) |
43 | { |
44 | register int i, j; |
45 | double *b, temp; |
46 | |
47 | /* Decompose matrix into L and U triangular matrices */ |
48 | if (lu_decompose(A, n) == 0) |
49 | return (0); /* Singular */ |
50 | |
51 | /* Invert matrix by solving n simultaneous equations n times */ |
52 | b = N_NEW(n, double); |
53 | for (i = 0; i < n; i++) { |
54 | for (j = 0; j < n; j++) |
55 | b[j] = 0.0; |
56 | b[i] = 1.0; |
57 | lu_solve(Ainv[i], b, n); /* Into a row of Ainv: fix later */ |
58 | } |
59 | free(b); |
60 | |
61 | /* Transpose matrix */ |
62 | for (i = 0; i < n; i++) { |
63 | for (j = 0; j < i; j++) { |
64 | temp = Ainv[i][j]; |
65 | Ainv[i][j] = Ainv[j][i]; |
66 | Ainv[j][i] = temp; |
67 | } |
68 | } |
69 | |
70 | return (1); |
71 | } |
72 | |