GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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>
9void 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 ldumat(double *a, double *u, int m, int n)
Definition ldumat.c:9