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