GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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"
10 int 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 }
struct ps_state ps
void free(void *)
double t
Definition: r_raster.c:39
int solv(double *a, double *b, int n)
Definition: solv.c:10
double b
Definition: r_raster.c:39