GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
minv.c
Go to the documentation of this file.
1/* minv.c CCMATH mathematics library source code.
2 *
3 * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
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"
10int 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}
double t
int minv(double *a, int n)
Definition minv.c:10
struct ps_state ps