GRASS Programmer's Manual  6.5.svn(2014)-r66266
minv.c
Go to the documentation of this file.
1 /* minv.c CCMATH mathematics library source code.
2  *
4  * This code may be redistributed under the terms of the GNU library
5  * public license (LGPL). ( See the lgpl.license file for details.)
6  * ------------------------------------------------------------------------
7  */
8 #include <stdlib.h>
9 #include "ccmath.h"
10 int minv(double *a, int n)
11 {
12  int lc, *le;
13
14  double s, t, tq = 0., zr = 1.e-15;
15
16  double *pa, *pd, *ps, *p, *q, *q0;
17
18  int i, j, k, m;
19
20  le = (int *)malloc(n * sizeof(int));
21  q0 = (double *)malloc(n * sizeof(double));
22  for (j = 0, pa = pd = a; j < n; ++j, ++pa, pd += n + 1) {
23  if (j > 0) {
24  for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
25  *q++ = *p;
26  for (i = 1; i < n; ++i) {
27  lc = i < j ? i : j;
28  for (k = 0, p = pa + i * n - j, q = q0, t = 0.; k < lc; ++k)
29  t += *p++ * *q++;
30  q0[i] -= t;
31  }
32  for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
33  *p = *q++;
34  }
35  s = fabs(*pd);
36  lc = j;
37  for (k = j + 1, ps = pd; k < n; ++k) {
38  if ((t = fabs(*(ps += n))) > s) {
39  s = t;
40  lc = k;
41  }
42  }
43  tq = tq > s ? tq : s;
44  if (s < zr * tq) {
45  free(le - j);
46  free(q0);
47  return -1;
48  }
49  *le++ = lc;
50  if (lc != j) {
51  for (k = 0, p = a + n * j, q = a + n * lc; k < n; ++k) {
52  t = *p;
53  *p++ = *q;
54  *q++ = t;
55  }
56  }
57  for (k = j + 1, ps = pd, t = 1. / *pd; k < n; ++k)
58  *(ps += n) *= t;
59  *pd = t;
60  }
61  for (j = 1, pd = ps = a; j < n; ++j) {
62  for (k = 0, pd += n + 1, q = ++ps; k < j; ++k, q += n)
63  *q *= *pd;
64  }
65  for (j = 1, pa = a; j < n; ++j) {
66  ++pa;
67  for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
68  *q++ = *p;
69  for (k = 0; k < j; ++k) {
70  t = 0.;
71  for (i = k, p = pa + k * n + k - j, q = q0 + k; i < j; ++i)
72  t -= *p++ * *q++;
73  q0[k] = t;
74  }
75  for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
76  *p = *q++;
77  }
78  for (j = n - 2, pd = pa = a + n * n - 1; j >= 0; --j) {
79  --pa;
80  pd -= n + 1;
81  for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
82  *q++ = *p;
83  for (k = n - 1, ps = pa; k > j; --k, ps -= n) {
84  t = -(*ps);
85  for (i = j + 1, p = ps, q = q0; i < k; ++i)
86  t -= *++p * *q++;
87  q0[--m] = t;
88  }
89  for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
90  *p = *q++;
91  }
92  for (k = 0, pa = a; k < n - 1; ++k, ++pa) {
93  for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
94  *q++ = *p;
95  for (j = 0, ps = a; j < n; ++j, ps += n) {
96  if (j > k) {
97  t = 0.;
98  p = ps + j;
99  i = j;
100  }
101  else {
102  t = q0[j];
103  p = ps + k + 1;
104  i = k + 1;
105  }
106  for (; i < n;)
107  t += *p++ * q0[i++];
108  q0[j] = t;
109  }
110  for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
111  *p = *q++;
112  }
113  for (j = n - 2, le--; j >= 0; --j) {
114  for (k = 0, p = a + j, q = a + *(--le); k < n; ++k, p += n, q += n) {
115  t = *p;
116  *p = *q;
117  *q = t;
118  }
119  }
120  free(le);
121  free(q0);
122  return 0;
123 }
tuple q
Definition: forms.py:2019
int minv(double *a, int n)
Definition: minv.c:10
void * malloc(YYSIZE_T)
void free(void *)
int n