GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-7d4cd0c030
csolv.c
Go to the documentation of this file.
1 /* csolv.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
9 #include <stdlib.h>
10 #include "ccmath.h"
11
12 int csolv(Cpx *a, Cpx *b, int n)
13 {
14  int i, j, k, lc;
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  q0 = (Cpx *)calloc(n, sizeof(Cpx));
23  pa = a;
24  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(q0);
54  return -1;
55  }
56  if (lc != j) {
57  h = b[j];
58  b[j] = b[lc];
59  b[lc] = h;
60  p = a + n * j;
61  q = a + n * lc;
62  for (k = 0; k < n; ++k) {
63  h = *p;
64  *p++ = *q;
65  *q++ = h;
66  }
67  }
68  t = pd->re * pd->re + pd->im * pd->im;
69  h.re = pd->re / t;
70  h.im = -(pd->im) / t;
71  for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
72  z.re = ps->re * h.re - ps->im * h.im;
73  z.im = ps->im * h.re + ps->re * h.im;
74  *ps = z;
75  }
76  }
77  for (j = 1, ps = b + 1; j < n; ++j, ++ps) {
78  for (k = 0, p = a + n * j, q = b, z.re = z.im = 0.; k < j; ++k) {
79  z.re += p->re * q->re - p->im * q->im;
80  z.im += p->im * q->re + p->re * q->im;
81  ++p;
82  ++q;
83  }
84  ps->re -= z.re;
85  ps->im -= z.im;
86  }
87  for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
88  for (k = j + 1, p = pd + 1, q = b + j + 1, z.re = z.im = 0.; k < n;
89  ++k) {
90  z.re += p->re * q->re - p->im * q->im;
91  z.im += p->im * q->re + p->re * q->im;
92  ++p;
93  ++q;
94  }
95  h.re = ps->re - z.re;
96  h.im = ps->im - z.im;
97  t = pd->re * pd->re + pd->im * pd->im;
98  ps->re = (h.re * pd->re + h.im * pd->im) / t;
99  ps->im = (h.im * pd->re - h.re * pd->im) / t;
100  --ps;
101  }
102  free(q0);
103  return 0;
104 }
int csolv(Cpx *a, Cpx *b, int n)
Definition: csolv.c:12
struct ps_state ps
double b
Definition: r_raster.c:39
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