Graphviz  2.41.20170921.2350
matinv.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
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 }
#define N_NEW(n, t)
Definition: memory.h:36
void lu_solve(double *x, double *b, int n)
Definition: lu.c:142
int matinv(double **A, double **Ainv, int n)
Definition: matinv.c:42
int lu_decompose(double **a, int n)
Definition: lu.c:67