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