GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
ldumat.c
Go to the documentation of this file.
1 /* ldumat.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 void ldumat(double *a, double *u, int m, int n)
10 {
11  double *p0, *q0, *p, *q, *w;
12 
13  int i, j, k, mm;
14 
15  double s, h;
16 
17  w = (double *)calloc(m, sizeof(double));
18  for (i = 0, mm = m * m, q = u; i < mm; ++i)
19  *q++ = 0.;
20  p0 = a + n * n - 1;
21  q0 = u + m * m - 1;
22  mm = m - n;
23  i = n - 1;
24  for (j = 0; j < mm; ++j, q0 -= m + 1)
25  *q0 = 1.;
26  if (mm == 0) {
27  p0 -= n + 1;
28  *q0 = 1.;
29  q0 -= m + 1;
30  --i;
31  ++mm;
32  }
33  for (; i >= 0; --i, ++mm, p0 -= n + 1, q0 -= m + 1) {
34  if (*p0 != 0.) {
35  for (j = 0, p = p0 + n, h = 1.; j < mm; p += n)
36  w[j++] = *p;
37  h = *p0;
38  *q0 = 1. - h;
39  for (j = 0, q = q0 + m; j < mm; q += m)
40  *q = -h * w[j++];
41  for (k = i + 1, q = q0 + 1; k < m; ++k) {
42  for (j = 0, p = q + m, s = 0.; j < mm; p += m)
43  s += w[j++] * *p;
44  s *= h;
45  for (j = 0, p = q + m; j < mm; p += m)
46  *p -= s * w[j++];
47  *q++ = -s;
48  }
49  }
50  else {
51  *q0 = 1.;
52  for (j = 0, p = q0 + 1, q = q0 + m; j < mm; ++j, q += m)
53  *q = *p++ = 0.;
54  }
55  }
56  free(w);
57 }
void free(void *)
void ldumat(double *a, double *u, int m, int n)
Definition: ldumat.c:9