GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
csolv.c
Go to the documentation of this file.
1/* csolv.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 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}
int csolv(Cpx *a, Cpx *b, int n)
Definition: csolv.c:10
double b
double t
struct ps_state ps
Definition: ccmath.h:38
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38