GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
cminv.c
Go to the documentation of this file.
1/* cminv.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 cminv(Cpx * a, int n)
11{
12 int i, j, k, m, lc, *le;
13
14 Cpx *ps, *p, *q, *pa, *pd;
15
16 Cpx z, h, *q0;
17
18 double s, t, tq = 0., zr = 1.e-15;
19
20 le = (int *)calloc(n, sizeof(int));
21 q0 = (Cpx *) calloc(n, sizeof(Cpx));
22 pa = pd = a;
23 for (j = 0; j < n; ++j, ++pa, pd += n + 1) {
24 if (j > 0) {
25 for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
26 *q++ = *p;
27 for (i = 1; i < n; ++i) {
28 lc = i < j ? i : j;
29 z.re = z.im = 0.;
30 for (k = 0, p = pa + i * n - j, q = q0; k < lc; ++k, ++q, ++p) {
31 z.re += p->re * q->re - p->im * q->im;
32 z.im += p->im * q->re + p->re * q->im;
33 }
34 q0[i].re -= z.re;
35 q0[i].im -= z.im;
36 }
37 for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
38 *p = *q++;
39 }
40 s = fabs(pd->re) + fabs(pd->im);
41 lc = j;
42 for (k = j + 1, ps = pd; k < n; ++k) {
43 ps += n;
44 if ((t = fabs(ps->re) + fabs(ps->im)) > s) {
45 s = t;
46 lc = k;
47 }
48 }
49 tq = tq > s ? tq : s;
50 if (s < zr * tq) {
51 free(le - j);
52 free(q0);
53 return -1;
54 }
55 *le++ = lc;
56 if (lc != j) {
57 p = a + n * j;
58 q = a + n * lc;
59 for (k = 0; k < n; ++k, ++p, ++q) {
60 h = *p;
61 *p = *q;
62 *q = h;
63 }
64 }
65 t = pd->re * pd->re + pd->im * pd->im;
66 h.re = pd->re / t;
67 h.im = -(pd->im) / t;
68 for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
69 z.re = ps->re * h.re - ps->im * h.im;
70 z.im = ps->im * h.re + ps->re * h.im;
71 *ps = z;
72 }
73 *pd = h;
74 }
75 for (j = 1, pd = ps = a; j < n; ++j) {
76 for (k = 0, pd += n + 1, q = ++ps; k < j; ++k, q += n) {
77 z.re = q->re * pd->re - q->im * pd->im;
78 z.im = q->im * pd->re + q->re * pd->im;
79 *q = z;
80 }
81 }
82 for (j = 1, pa = a; j < n; ++j) {
83 ++pa;
84 for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
85 *q++ = *p;
86 for (k = 0; k < j; ++k) {
87 h.re = h.im = 0.;
88 for (i = k, p = pa + k * n + k - j, q = q0 + k; i < j; ++i) {
89 h.re -= p->re * q->re - p->im * q->im;
90 h.im -= p->im * q->re + p->re * q->im;
91 ++p;
92 ++q;
93 }
94 q0[k] = h;
95 }
96 for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
97 *p = *q++;
98 }
99 for (j = n - 2, pd = pa = a + n * n - 1; j >= 0; --j) {
100 --pa;
101 pd -= n + 1;
102 for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
103 *q++ = *p;
104 for (k = n - 1, ps = pa; k > j; --k, ps -= n) {
105 z.re = -ps->re;
106 z.im = -ps->im;
107 for (i = j + 1, p = ps + 1, q = q0; i < k; ++i, ++p, ++q) {
108 z.re -= p->re * q->re - p->im * q->im;
109 z.im -= p->im * q->re + p->re * q->im;
110 }
111 q0[--m] = z;
112 }
113 for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
114 *p = *q++;
115 }
116 for (k = 0, pa = a; k < n - 1; ++k, ++pa) {
117 for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
118 *q++ = *p;
119 for (j = 0, ps = a; j < n; ++j, ps += n) {
120 if (j > k) {
121 h.re = h.im = 0.;
122 p = ps + j;
123 i = j;
124 }
125 else {
126 h = q0[j];
127 p = ps + k + 1;
128 i = k + 1;
129 }
130 for (; i < n; ++i, ++p) {
131 h.re += p->re * q0[i].re - p->im * q0[i].im;
132 h.im += p->im * q0[i].re + p->re * q0[i].im;
133 }
134 q0[j] = h;
135 }
136 for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
137 *p = *q++;
138 }
139 for (j = n - 2, le--; j >= 0; --j) {
140 for (k = 0, p = a + j, q = a + *(--le); k < n; ++k, p += n, q += n) {
141 h = *p;
142 *p = *q;
143 *q = h;
144 }
145 }
146 free(le);
147 free(q0);
148 return 0;
149}
int cminv(Cpx *a, int n)
Definition: cminv.c:10
double t
struct ps_state ps
Definition: ccmath.h:38
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38