GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
9#include <stdlib.h>
10#include "ccmath.h"
11
12int 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}
struct complex Cpx
Definition ccmath.h:41
int csolv(Cpx *a, Cpx *b, int n)
Definition csolv.c:12
double b
double t
struct ps_state ps
double re
Definition ccmath.h:39
double im
Definition ccmath.h:39