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"
39extern int lu_decompose(double **a, int n);
40extern void lu_solve(double *x, double *b, int n);
41
42int 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