GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
house.c
Go to the documentation of this file.
1/* house.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#include "ccmath.h"
10void house(double *a, double *d, double *dp, int n)
11{
12 double sc, x, y, h;
13
14 int i, j, k, m, e;
15
16 double *qw, *qs, *pc, *p;
17
18 qs = (double *)calloc(2 * n, sizeof(double));
19 for (j = 0, qw = qs + n, pc = a; j < n; pc += n + 1)
20 qw[j++] = *pc;
21 for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
22 m = n - j - 1;
23 for (i = 1, sc = 0.; i <= m; ++i)
24 sc += pc[i] * pc[i];
25 if (sc > 0.) {
26 sc = sqrt(sc);
27 if ((x = *(pc + 1)) < 0.) {
28 y = x - sc;
29 h = 1. / sqrt(-2. * sc * y);
30 }
31 else {
32 y = x + sc;
33 h = 1. / sqrt(2. * sc * y);
34 sc = -sc;
35 }
36 for (i = 0, qw = pc + 1; i < m; ++i) {
37 qs[i] = 0.;
38 if (i)
39 qw[i] *= h;
40 else
41 qw[i] = y * h;
42 }
43 for (i = 0, e = j + 2, p = pc + n + 1, h = 0.; i < m;
44 ++i, p += e++) {
45 qs[i] += (y = qw[i]) * *p++;
46 for (k = i + 1; k < m; ++k) {
47 qs[i] += qw[k] * *p;
48 qs[k] += y * *p++;
49 }
50 h += y * qs[i];
51 }
52 for (i = 0; i < m; ++i) {
53 qs[i] -= h * qw[i];
54 qs[i] += qs[i];
55 }
56 for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
57 for (k = i; k < m; ++k)
58 *p++ -= qw[i] * qs[k] + qs[i] * qw[k];
59 }
60 }
61 d[j] = *pc;
62 dp[j] = sc;
63 }
64 d[j] = *pc;
65 dp[j] = *(pc + 1);
66 d[j + 1] = *(pc + n + 1);
67 for (j = 0, pc = a, qw = qs + n; j < n; ++j, pc += n + 1) {
68 *pc = qw[j];
69 for (i = 1, p = pc + n; i < n - j; p += n)
70 pc[i++] = *p;
71 }
72 free(qs);
73}
void house(double *a, double *d, double *dp, int n)
Definition house.c:10
#define x