GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
solv.c
Go to the documentation of this file.
1/* solv.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 general
5 * public license. ( See the gpl.license file for details.)
6 * ------------------------------------------------------------------------
7 */
8#include <stdlib.h>
9#include "ccmath.h"
10int solv(double *a, double *b, int n)
11{
12 int i, j, k, lc;
13
14 double *ps, *p, *q, *pa, *pd;
15
16 double *q0, s, t, tq = 0., zr = 1.e-15;
17
18 q0 = (double *)calloc(n, sizeof(double));
19 for (j = 0, pa = a, pd = a; j < n; ++j, ++pa, pd += n + 1) {
20 if (j) {
21 for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
22 *q++ = *p;
23 for (i = 1; i < n; ++i) {
24 lc = i < j ? i : j;
25 for (k = 0, p = pa + i * n - j, q = q0, t = 0.; k < lc; ++k)
26 t += *p++ * *q++;
27 q0[i] -= t;
28 }
29 for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
30 *p = *q++;
31 }
32 s = fabs(*pd);
33 lc = j;
34 for (k = j + 1, ps = pd; k < n; ++k) {
35 if ((t = fabs(*(ps += n))) > s) {
36 s = t;
37 lc = k;
38 }
39 }
40 tq = tq > s ? tq : s;
41 if (s < zr * tq) {
42 free(q0);
43 return -1;
44 }
45 if (lc != j) {
46 t = b[j];
47 b[j] = b[lc];
48 b[lc] = t;
49 for (k = 0, p = a + n * j, q = a + n * lc; k < n; ++k) {
50 t = *p;
51 *p++ = *q;
52 *q++ = t;
53 }
54 }
55 for (k = j + 1, ps = pd, t = 1. / *pd; k < n; ++k)
56 *(ps += n) *= t;
57 }
58 for (j = 1, ps = b + 1; j < n; ++j) {
59 for (k = 0, p = a + n * j, q = b, t = 0.; k < j; ++k)
60 t += *p++ * *q++;
61 *ps++ -= t;
62 }
63 for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
64 for (k = j + 1, p = pd, q = b + j, t = 0.; k < n; ++k)
65 t += *++p * *++q;
66 *ps -= t;
67 *ps-- /= *pd;
68 }
69 free(q0);
70 return 0;
71}
double b
double t
struct ps_state ps
int solv(double *a, double *b, int n)
Definition solv.c:10